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Abstract 

In a body periodically strained by tides, heating produced by viscous friction is 
far from homogeneous. The spatial distribution of tidal heating depends in a com- 
plicated way on the tidal potential and on the internal structure of the body. I show 
here that the distribution of the dissipated power within a spherically stratified body 
is a linear combination of three angular functions. These angular functions depend 
only on the tidal potential whereas the radial weights are specified by the internal 
structure of the body. The 3D problem of predicting spatial patterns of dissipation 
at all radii is thus reduced to the ID problem of computing weight functions. I com- 
pute spatial patterns in various toy models without assuming a specific rheology: a 
viscoelastic thin shell stratified in conductive and convective layers, an incompress- 
ible homogeneous body and a two-layer model of uniform density with a liquid or 
rigid core. For a body in synchronous rotation undergoing eccentricity tides, dissi- 
pation in a mantle surrounding a liquid core is highest at the poles. Within a soft 
layer (or asthenosphcrc) in contact with a more rigid layer, the same tides generate 
maximum heating in the equatorial region with a significant degree-four structure if 
the soft layer is thin. The asthenosphere can be a layer of partial melting in the up- 
per mantle or, very differently, an icy layer in contact with a silicate mantle or solid 
core. Tidal heating patterns are thus of three main types: mantle dissipation (with 
the icy shell above an ocean as a particular case), dissipation in a thin soft layer and 
dissipation in a thick soft layer. Finally, I show that the toy models predict well pat- 
terns of dissipation in Europa, Titan and lo. The formalism described in this paper 
applies to dissipation within solid layers of planets and satellites for which internal 
spherical symmetry and viscoelastic linear rheology are good approximations. 
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1 Introduction 



Fifty years ago, William Kaula found that the heat dissipated by tidal friction within 
the Moon is extremely nonuniform both radially and laterally \Kaula 1963 1964 . At 



that time, nonuniform tidal heating was a rather academic subject as these variations 
were not observable: tidal heating in the Moon is indeed much smaller than radiogenic 
heating. Spatial variations of tidal heating were actually a byproduct of computing the 
total power dissipated in the body, a key factor in modeling orbital evolution. Kaula's 
calculations were based on the microscopic (or micro) approach to tidal dissipation, 
which starts with the computation of viscoelastic tidal strains. A bit earlier, Munk an^ 
MacDonald 19601 had given an approximate formula for the total power dissipated by 



tides with a macroscopic (or macro) approach. This alternative approach does not re- 
quire the computation of tidal strains: the tidal bulge lags the tidal forcing by an angle 
parameterizing the viscous response. Their formula, however, was limited to an incom- 
pressible homogeneous body, in contrast with the micro approach which is applicable to 
any model of internal structure. 

New theoretical developments had to wait until 1978, when two papers significantly 



improved tidal heating computations. First, Peale and Cassen 11978 reconsidered both 



macro and micro approaches to tidal heating, correcting several errors in the various 
formulas. They also mapped tidal heating variations for eccentricity tides and obliquity 
tides, showing that tidal dissipation in a homogeneous Moon is maximum at the poles 
in the former case and at the equator in the latter. Furthermore they showed that the 
presence of a large liquid core enhances dissipation in the mantle. Second, Zschau [1978 



derived a simple formula for the total dissipated power in a spherically stratified com- 
pressible body, in which the influence of the body's internal structure appears through 
Im{k2), the imaginary part of the gravity tidal Love number k2- Zschau's formula was 
the first step toward reconciling the macro and micro approaches, since k2 can be com- 



puted if the internal structure of the body is known. Tobie et al. 2005b| later applied 
the variational method to relate /m(/c2) to the imaginary part of the volume- integrated 
strain power. They also computed the radial distribution of the dissipated power in 
terms of deformation functions that can be evaluated with standard methods for any 



spherically stratified model. I will show that the formulas of Zschau 1978 and Tobie 



et al. 2005b are recovered by spatially averaging the local power obtained in the micro 
approach, thus bridging the last gap between the macro and micro approaches. My 
paper, however, is primarily about the angular distribution of the dissipated power. 

Spatial variations in tidal heating became relevant when spacecraft sent back incred- 
ible surface data showing that the Galilean satellites lo and Europa undergo strong tidal 
deformations and heating. Tidal heating was the only explanation for lo's volcanism 



Peale and Cassen, 1978 but what was going on beneath the surface was a mystery. 



Maybe the distribution of volcanoes could be used in order to constrain the internal 
structure of the satellite? Segatz et al. 1988] suggested that dissipation occurred ei- 



ther in the whole mantle or mostly in a thin asthenosphere close to the surface. These 
two models predict completely different patterns of surface heat fiux with maximum 
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dissipation at the poles in the former case and at the equator in the latter (a mix of 
the two is of com'se possible). Galileo data have often been interpreted as favoring the 



asthenospheric dissipation 


model though the issue remains controversial Lopes- Gautier 


et al. , 1999 ; Tackley et al. 


2001 


Kirchoffet al. 2011 


Veeder et al. , 2012 Hamilton et al. 


2012 . Another idea consists in using long wavelength topography as an indirect measure 



of the heat flux (valid if the topography is isostatically compensated) , but the promising 



analysis of Voyager data Ross et al. , 1990 was not confirmed by Galileo measurements 
\Thomas et al. 1998 . Lateral variations of surface heat flux however also depend on the 
heat transport mechanism which could be determinant. 

On icy satellites like Europa, Titan and Enceladus, spatial variations of tidal heating 
are interesting for other reasons. In these satellites, tidal heating can melt the ice at 
depth and create a global subsurface ocean. The covering icy shell varies in thickness be- 
cause of spatial variations in tidal heating and solar insolation. Ojakangas and Stevenson\ 
1989a|b computed icy shell thickness variations on Europa with the aim of predicting 



Nimmo et al. 20071 and Nimmo and Bills 



nonsynchronous rotation and polar wander 
[2010 used the same method to predict long wavelength topography on Europa and on 
Titan, respectively. Variations of the thickness of the icy shell (or more generally the 



lithosphere) also influence surface tectonics \Beuthe 2010 



Finally, spatial variations of tidal heating are important because of their strong 
coupling to convection. Several convection models have used as input tidal heating 
predicted by spherically stratified models | Tackley et al. 2001 ; Tobie et al. , 2003 Roberts 



and Nimmo , 2008 . An important limitation of this approach is the neglect of lateral 
viscosity variations on tidal heat production, which can be taken into account by giving 
up the assumption of spherical symmetry and solving simultaneously for convection and 
tidal dissipation Behounkovd et al.. 2010; Han and Showman. 2010 . Chaos terrain on 



Europa could be a visible result of spatially varying tidal heating enhanced by a local 
drop in viscosity. 

Until now, predicting spatial patterns of tidal dissipation meant computing the dis- 
sipated power at every point within the body. This laborious procedure obscures the 
link between the internal structure and the resulting pattern, and makes it difficult to 
look for all possible patterns generated by realistic internal structures. In this paper, 
I show that the dissipated power at a given radius within a spherically stratified body 
is the linear combination of three basic patterns that depend only on the tidal poten- 
tial. The coefficients weighting the patterns depend are radial functions which can be 
computed with standard methods developed for tidal deformation problems once the 
internal structure of the body has been specified. 

I study the infiuence of the internal structure on dissipation patterns by computing 
the dissipation weight functions in various toy models without assuming a specific rhe- 
ology. The toy models are the thin icy shell above an ocean, the homogeneous body 
(relevant to a completely solid body with little stratification or to a solid core) and the 
incompressible two- layer body, either with a liquid core or with rigid core, the latter 
case being relevant to dissipation in a soft layer (such as an asthenosphere) above a 
more rigid layer. It is well-known that dissipation patterns are completely different if 
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dissipation occurs in the deep mantle and in a thin asthenosphere. Besides these two 
classes of patterns, I show that dissipation in a thick asthenosphere leads to a third 
type of dissipation pattern, with maximum heating at the equator as in asthenospheric 
dissipation but with a lower content in harmonic degree four. Toy models predict well 
dissipation patterns within real bodies though not their magnitude. I will give three ex- 
amples of this by computing dissipation weight functions for realistic internal structures 
of Europa, Titan and lo. 



2 Dissipated power 

2.1 Power, strains and tidal potential 

In the micro approach to dissipation, the dissipated power within the planet or satellite is 
expressed in terms of tidal strains (this formula is derived in Appendix A and compared 
with other expressions found in the literature). If tides operate at only one angular 
frequency u, the dissipated power per unit volume averaged over one orbital period is 
given by 

3"^" 



P 



3lm{ji) 



+ '^Im{K)\ef 



(1) 



where jl (resp. K) is the complex shear (resp. bulk) modulus, e^^ is the Fourier transform 
of the strain tensor and e is the trace of e^j. All quantities implicitly depend on the 
frequency lv and on the point x within the planet where the power is evaluated. In 
this paper, the tilde on viscoelastic parameters indicates that they are complex and 
frequency-dependent (the tilde is dropped if the parameters are purely elastic). 

If there are several tidal frequencies, the total power averaged over time is a sum 



over these frequencies (interferences vanish, see Eq. (68)) so that it becomes essential to 
know the frequency dependence of the viscoelastic parameters (i.e. the rheology). Unfor- 
tunately the rheology of planetary bodies is poorly constrained \Jackson 2007 Karato 



2008 . Earth's mantle has been mainly studied at frequencies that are much higher (lab- 



oratory experiments), moderately higher (seismic attenuation and seismic anisotropy) or 



much lower (Chandler wobble and postglacial rebound) than tidal frequencies \Karato 



2010 


|. It is indeed difficult to determine Earth's viscous 


[e.g. 


Benjamin et al. 


2006 




Nakada and Karato 




2012 and 



Maxwell rheology is often used in studies of tidal dissipation because it is the simplest 
model in which the response changes from elastic to viscous as the frequency decreases. 
It is however not clear how the Maxwell viscosity is related to the true viscosity of the 



material Ross and Schubert 


1986 


Bills et al. 




2005 




Sotin et al. 




2009 . Rheological 


models depending on more parameters such as the Andrade model | 


Castillo-Rogez et al. 



2011 or the extended Burgers model iNimmo et al. , 2012 could be more realistic. More- 



over there has been a long-standing debate on whether viscous deformations in Earth's 
mantle are mainly due to diffusion creep or to dislocation creep, corresponding to a lin- 
ear or nonlinear rheology, respectively Karato and Wu, 1993 . In this paper, I assume 



that the rheology is linear without being more specific about it except in applications 
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to real bodies for which I use the Maxwell model. Besides I consider only the dominant 
tidal frequency; this restriction is appropriate for Solar System bodies, but it should be 
lifted for exoplanets with short orbital periods. 

Dissipation within the Earth is generally much larger in shear than in uniform com- 

. For 



pression {Im{jl) ^> Im{K)), at least in the seismic frequency band Anderson 



1989 



this reason, most models of planetary dissipation assume that the bulk modulus is real 



and independent of frequency \Ranalli 1995, p. 142], though it is not necessarily true in 
presence of a high fraction of partial melt Schmeling , 1985 . In Earth's asthenosphere. 



the ratio of bulk to shear dissipation could possibly reach 30% Durek and Ekstrom 



1995 



Bulk dissipation should thus be kept in mind when studying bodies (such as lo) 
where asthenospheric dissipation could be the dominant mechanism. 

The strains appearing in Eq. ([T]) are induced by tidal deformations which are in turn 
related to the tidal potential. The tidal potential at the surface of the deformed body 
can be expanded in spherical harmonics of degree £ and order m \Kaula 1964 , each 



component being the superposition of several terms of angular frequencies w^^j. 



U{t,9,<f>) = Re 



where Re{x) means real part of x, 9 is the colatitude and (j) is the longitude in a frame 
attached to the body. The coefficients ^imj{di4>) are defined by this equation if you 
know the tidal potential from the formulas in Kaula [1964| ; they are not the same as 
the complex coefficients of the Fourier series of U{t, 9, 4>) since the frequencies uimj are 
not all different. If the orbital eccentricity and the obliquity of the body are not zero, 
an infinite number of terms (indexed by j) contribute to the potential at a given degree 
i and order m. I neglect all tidal perturbations except the dominant contribution of 
degree i = 2 and I sum over the order m. The component of the tidal potential at 
frequency a; can be written as 



U^{t,9,(l)) 



i?e($ uj{9,(l)) exp (iwt)) 



exp{iujt) + c.c. - 



(3) 



where c.c. means complex conjugate. 

Strains are related to the tidal potential through displacements. If the internal 
structure of the body is spherically symmetric, displacements due to tides can be written 
as 



Alterman et al. , 1959 , Takeuchi and Saito , 1972 Saito , 1974 



Ur 


= Vi 




Ug 


= Vs 


^"■^ d9 ' 




= 2/3 


. . 1 d'^^ 

''^^ sm9 86 



(4) 



where {ur, ug, u^) are the Fourier components of the displacement defined as in Eqs. (62) 
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The two radial functions yi(r) and ysir) depend on the internal structure of the 
body and have the dimension of inverse acceleration. These functions are part of a set 
of six radial functions {yi, ...^yo) that solve six coupled linear differential equations of 
first order (I drop from now on the explicit dependence on r). Strains depend on yi, 
2/3 and their derivatives. If derivatives are eliminated with the help of the differential 
equations, strains also depend on the functions y2 and y^ respectively associated with 
the stresses arr and cxre (or cr-r^)- The remaining functions ys and yg are related to 
the gravity potential. The derivatives {y'i,y'^) are related to {yi,...,y^) by Eqs. (82) of 
Takeuchi and Saito [1972 : 

1 ^ry2 -{K- 2/i/3) (2yi - Gyg) 



ry'i 



K + 4/2/3 
ry'z - 2/3 + yi • 



The quantity y^^/fi actually measures the shear strain (see Eqs. 
pressible medium, the first equation reduces to 



ry'i 



-2yi + 6y3 



Since stresses with a radial component vanish at the surface r 
conditions for y2 and y^ are 

y^{R)=y^{R)=Q. 



(5) 

^ below). In an incom- 
(6) 

R, the boundary 
(7) 



1972 


)■ I 


2004 


are 


1972 


• 


Kaula 



Beware that the definitions of yi vary between authors (I follow here Takeuchi and Saito 
articular, the functions (yi, y2, ys, y45 ys, ye) oi \Sahadini and Vermeersen 
Liivalent to the functions (yi, ys, y2, y4, — ys, — ye) of Takeuchi and Saito 



Kaula [1963 1964 already used the y^ functions in his pioneering papers on tidal 



dissipation; his example was recently followed by Tabic et al. |2005b and Roberts and 



Nimmo [2008 . This formalism has the advantage that it is widely used in geophysics to 



compute the deformation of Earth's surface for complicated internal structures \Spada 



et al. 



1911 



2011 . Peale and Casscn 1978 



unfortunately returned to the method of Love 
which is difficult to apply correctly to structures more complicated than two 
layers of the same density. 

In spherical coordinates, I divide strain components in three classes: (1) the purely 
radial component e^r is the radial strain, (2) the purely angular components {cqq, e^p^p, eg^) 
are the tangential strains, and (3) the mixed components (e^-e, Crtp) are the radial-tangential 
shear strains. I now insert Eqs. (|4|) into the strain-displacement equations in spherical 
coordinates 



Takeuchi and Saito 




1972 


Eqs. (20)]. Note that the non-diagonal strains 


in 


Takeuchi and Saito 


1972 must be multiplied by 1/2 if they are to be 



the components of the strain tensor. The radial and shear strains read 

err = y'l ^LU , 

ly4 d^gj 
2 ft de ' 
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2 fl sine dcp ' ^ ^ 



Doing the same for the tangential strains, I get 



= \ (^3 ^2 + yi) , 
ee<j> = -ygOs^a;, (9) 



where Oi 2 3 are differential operators of degree two on the sphere defined by Eqs. (80)- 



(81). 



Though Eqs. (|8])-([9]) are well-known, their tensorial structure with respect to rota- 
tions is always ignored in geophysics: e^r is a scalar, Ira are the components of a vector 
and lap are the components of a second-order symmetric tensor (Greek indices denote 
6 01 (j)). The tensorial character of strains under rotations is the key to simplifying the 
formula for the dissipated power. 

2.2 Strain invariants 

From a technical point of view, this is the key section of the paper. In order to evaluate 
the dissipated power in terms of the tidal potential, I need to compute in Eq. ([T]) the 
quantities |e|^ and e^j e*j which are invariant under arbitrary coordinate transformations 



Malvern 1969 . This general invariance is not obvious once the strains have been ex- 
pressed in spherical coordinates. Nevertheless invariance under rotations should remain 
obvious in that basis. It is indeed easy to construct rotationally invariant terms: take 
the product of two scalars (such as err), take the scalar product of two vectors (such as 
ira), contract two tensors (such as iaf^), or take the trace of a tensor. 



Let us define 



IT' 2 |~|2 

Ediiat = r \e\ 



E2 = r^e^^.e*., (10) 

with the subscript dilat denoting dilatation^ since e is the infinitesimal change of volume. 
In spherical coordinates, Ediiat reads 

2 I ~ ~ ~ 1 2 

Ediiat =1' \^rr + ^ee + ^(t>(t>\ ■ (H) 

The terms irr (scalar) and ioe + e^^ (trace of a tensor) are each invariant under rotations 



but it is easier to keep them together in Eq. (11). E2 is the sum of three terms that are 
each invariant under rotations: 

E2 — Eradial ~l~ -^s/iear ~l~ Etangent i (1^) 
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where 



^radial 
■^shear 



2 ~ ~* 



E, 



tangent 



2r^ {e^Q e*0 + e^^ e*^) , 



(13) 



Eradiai, E shear and Etangent are invariant because they are respectively the product of 
two scalars, the scalar product of two vectors and the contraction of two tensors. 

Inserting Eqs. ([8])-([9]) and the first Eqs. (82) into Eq. (11), I write the first invariant 

as 

EdUat = I {ry[ + 2yi + A) |^ (14) 

where A is the spherical Laplacian defined by Eq. ( |76[ ). Inserting Eqs. ([s]) and (79) into 
the first and second Eqs. (13), I write the radial and shear invariants as 



Eradiai 
Eghear 



I ' |2 l/is |2 



(15) 



Inserting Eqs. ^ and (82) into the third Eqs. (13), I write the tangential invariant as 
the sum of two terms invariant under rotations: 



E, 



tangent 



E, 



tanl 



+ E, 



tan2 • 



(16) 



where 



Etanl 
Etan2 



2|(yi+y3A/2)<D^|2 . 



(17) 



Noting that <I>^ is a spherical harmonic degree two, I insert Eqs. (77) and (83) into 
Eqs. pll) to (fTTl): 



Edilat 
Eradiai 
Eshear 
Etanl 
Etan2 



2 

I ' |2 1^ |2 



2 



2 



(l/4)|rV/ir(A-252)|^, 
(1/4) lysl' (AA - 2 (2^2 + 1) A + 2,^2 (62 + 2)) l^^f 
2\y, + i62/2)yf\^^ 



2 



(18) 



where 62 = —6 is the degree-two eigenvalue of the spherical Laplacian A (see Eq. (|77|). 
Strain invariants are now factorized into radial and angular parts. 
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2.3 Factorized power 

The last step consists in rewriting the dissipated power given by Eq. ([T]) in terms of the 
various strain invariants of Section 12.21 



OJ 



Im{jl) Eradial + Eg hear + Etanl + E^ 



tan2 



EdUat ] + ^ Im{K) Ediiat . (19) 



Examining the factorized strain invariants (Eqs. (|18|)), we see that the power at a given 
radius r depends on the angles {9, (p) through a linear combination of three factors, 



showing that there are at most three independent spatial patterns. 

2 

The three invariants without derivatives of \^uj\ can be combined as 

Eradial + Etan2 - EdUat/^ = (2/3) - Vi - {62/2) ^3 \^uj\ 



(20) 



(21) 



The three functions ( 20 ) are not the best choice of basic angular functions because of 
compensations between the three terms. Instead I will maintain as much as possible the 
distinction between the radial, radial-tangential shear and tangential dissipative terms. 
These contributions have indeed very different magnitudes in very different interior mod- 
els (see Section 4.3). Besides |<I>a;|^5 I thus propose to use as basic angular functions the 
combinations appearing in E^hear ^ind Efani (see Eqs. (18)): 

|2 



ni?)~*|$^, 
^ (A + I2)vl/^, 



12 
1 

48 



(AA + 22A + 48) . 



(22) 



These functions are dimensionless since the factor {nR)~^ absorbs the dimension of the 
squared tidal potential (n is the mean motion and R is the surface radius, see Eq. ([89])). 
The numbers 12, 22 and 48 result from setting 82 = —6 in Eqs. (18). The patterns 
and ^ are normalized so that the terms without derivatives are the same in the three 
functions. In the following, the index J collectively refers to {A^B^C). 



Inserting Eqs. ( |18| ), (21) and (22) into Eq. (19), I write the dissipated power as 
uj{nRf 



P 



2r2 



/m(/i) (/^ + /b + fc ^c) + MK) Hj, 



(23) 



where the weight functions fj and Hj^ are positive and depend on the internal radial 
functions y^: 

2 

6 



fA 



Vy'i 



c 



H 



K 



24|y3l' , 



(24) 
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The weights functions have the dimension of squared inverse acceleration. When aver- 
aging over angles, I will need the sum of the weight functions fj denoted as 



H^ = fA + fB + fi 



c • 



(25) 



In practice, y[ is evaluated from y^, and (see Eqs. ([s])). The terms B, C and K 
are associated with radial-tangential shear, tangential and bulk dissipation, respectively, 
whereas the term ^ is a combination of radial and tangential dissipation. The terms A 
and K depend on the same angular function and can be grouped into one term. Formulas 



(22), (23) and (24) are the central results of this paper. 



At the surface, weight functions can be related to the displacement Love numbers 
which characterize how the body responds to tidal forcing. The displacements yi{R) 
and yz{R) are indeed proportional to the tidal Love numbers /12 and I2, while 2/2 (-R) and 
y4{R) vanish because of Eqs. ([T]): 



ho h 

— ,o,-,o 

9 9 



(26) 



where g is the surface gravity. Moreover, the first of Eqs. ^ evaluated at the surface 
gives 

D 2P /12-3/2 , . 

-Ryi(^) = -Y3^ — - — ' (27) 

where I replaced fl and K by Young's modulus E and Poisson's ratio D: 



K 



E 



3(1 - 2D) 
E 

2(1 + z>) ■ 



(28) 



I compute the weight functions at the surface with Eqs. (24) to (27): 



if A t Ib ■• fc . 



H 



K 



r=R 



r=R 




h2 



0,24 



h2 



ho 



(29) 



The amplitude of tidal deformation (quantified by I/12I) varies by orders of magnitude 
depending on the global structure of the body. By contrast, the ratio 1/2/^2! typically 



ranges from 0.2 to 0.3 \Beuthe 2010, Fig. 10]. Eqs. (29) show that ratios of weight 



functions at the surface are insensitive to the magnitude of the viscoelastic response 
|/i2|- In good approximation, this property also holds inside the body. It is thus useful 
to define dimensionless rescaled weight functions by 



ho 



(30) 
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For example, rescaled weight functions do not depend on the rheology if the body is 



incompressible and homogeneous (Section 4.2) or if it is consists of an incompressible 



viscoelastic mantle above a liquid core of the same density (Section 4.3). 



Weight functions have a few properties that do not depend on a specific internal struc- 



ture. First, always vanishes at the surface (see Eqs. (29)) and at internal fluid/solid 
interfaces (where a condition similar to Eq. ([7]) holds). Pattern B can thus be disre- 
garded close to such boundaries, for example in a thin icy shell above an ocean. Second, 
fj^ and vanish in the incompressible limit at the boundary of an infinitely rigid layer 
because Hi = = there. Dissipation thus approximately follows Pattern B close to 
rigid/ viscoelastic boundaries, for example at the bottom of an icy shell in contact with 



a silicate mantle. Third, I can relate Hj^ to at the surface with Eqs. (29): 



Hk 




4 


A 


2 

= 3 


1 - 2D 


fA 


r=R 


3 


K 




l + D 



(31) 



This ratio is less than one- half for silicates {v ~ 1/4) and less than one- fifth for ice 
(j^ ~ 1/3). If the material is incompressible, it is possible to go further by applying 
Eq. 



fA 



I2|yi 
0, 



3^3 



(32) 



so that Hj^ is expected to be much smaller than at all radii if incompressibility is a 
good approximation. Since dissipation also receives contributions from terms B and C, 
the contribution of bulk dissipation to total dissipation remains in general minor even if 
Im{K) Im{fl) (see discussion in Section^. 



3 Spatial patterns, global power and surface flux 
3.1 Harmonic content of angular functions 



The angular functions t defined by Eqs. (22) can be computed from the expansion of 



the squared norm of the tidal potential in spherical harmonics. Since the tidal potential 
is of degree two, the expansion of its squared norm only includes functions of degrees 
zero, two and four: 

|<I>^|2 = (ni?)4(M/o + ^2 + ^4) , (33) 

where 

e 

^i=J2 ^i^rn (cos 61) (o^^ COS m(j) + bg^ sin m(j)) . (34) 

m=0 

The constant = Oqq is the spatial average of the squared norm of the nondimensional 
tidal potential, 

= i^^r^)— f l^^fsinB dBdcl), (35) 
An , q 



3 SPATIAL PATTERNS, GLOBAL POWER AND SURFACE FLUX 



12 



Table 1: Squared norm of the nondimensional tidal potential for synchronous rotation: non- 
zero spherical harmonic coefficients (a 
obliquity tides (e 
s„ and c. 



j) of Eq. (34) for eccentricity tides (/ 
~ ~ (|93l. 



0). The coefficients are computed with Eqs. (89) and (92) 
denote sinwp and coswp, respectively, where LOp is the argument of pericentre 



= 0) and for 
The symbols 







'^20 


^^22 






044 


eccentricity tides (x e^) 


21/5 


-33/7 


9/14 


387/140 


-27/140 


-3/160 


obliquity tides (x s'm^I) 


3/5 


3/7 


3/14 


-36/35 


3/35 







^21 


«41 


"43 


^21 


&41 


^43 


interference ( x e sin /) 


(6/7)sp 


(-81/70)sp 


(9/140)sp 


(12/7)cp 


(-18/35)cp 


(3/35)cp 



while ^'2 and ^'4 represent spatial variations with zero average. Table [T] gives the nu- 
merical values of the coefficients for a satellite in synchronous rotation with non-zero 
eccentricity and obliquity. 

In terms of the harmonic functions ^'£, the angular functions read 

*B = ^0 + 1^2-1^4, 

= ^o-*2+g^4- (36) 
The angular dependence cancels in the following combination: 

^^ + 2^'5 + 2^c = 5^o- (37) 



Eqs. (36) give us a qualitative understanding of the spatial patterns. ^'^ and mainly 
differ by the sign of the degree-four term. By contrast, the term of degree two in is 
of opposite sign with respect to ^'^ (or ^'g) and the degree-four contribution to "if^j is 
much smaller, even though contains derivatives of the fourth order. The harmonic 



content of the angular functions can be measured with the variance (see Eq. (94)) if one 
knows the tidal potential. Table [2] shows that the content of the variance in degree four 
is much larger for obliquity tides than for eccentricity tides though is still mainly of 
degree two. Fig. [T] shows the basic patterns of tidal heating (angular functions ^'j) for 
a body in synchronous rotation with either eccentricity tides or obliquity tides, drawn 
from Eqs. (36) in which I substitute Eq. (34) and the values of Table [T] 



The relative importance of the three patterns can be determined from the comparison 
of the weight functions but one should keep in mind that the do not have the same 
variance. I thus define equal-variance weight functions fj through 

'^ = jB7)f^ (J = 4.B.C). (38) 

The standard deviations sd(^'j) are tabulated in Table |3] for synchronous rotation (ec- 
centricity tides only or obliquity tides only). For eccentricity tides, fg/f^ = 0.56 and 
fcllc = 0-84. For obliquity tides, fs/fj, = fc/fc = 0-61- 
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Table 2: Harmonic content of the basic spatial patterns for synchronous rotation: contribution 
of degree four (in %) to the variance of for eccentricity tides (/ = 0) and for obliquity tides 
(e = 0). The remainder of the variance is due to degree two. 









*c 


unit 


eccentricity tides 


30 


44 


1 


% 


obliquity tides 


64 


76 


5 


% 



Table 3: Standard deviation of the basic spatial patterns for synchronous rotation. 



eccentricity tides (x e^) 
obliquity tides (x sin^I) 



'f = 2.80 



= 0.64 



f = 1.55 



= 0.39 



2.35 



= 0.39 



Instead of \l'j, I could use as basic angular functions and express the spatially 
varying part of the power as a linear combination of only two functions (\I'2 and 'I'4). 
The radial weights associated with ^'2 and ^4, however, are not always positive. For 
example, the degree- two component of the power has a minimum either at the poles or 
at the equator, depending on the sign of the radial function. Physically, these two cases 
are interpreted as different spatial patterns. 



3.2 Global power 



I show here that the dissipated power integrated over the volume of the body is equal 
to the global power computed with the macro approach to tidal heating discussed in 
the Introduction. The angular average of the local power at radius r can be read from 
Eqs. ([23l), ([25l) and ([36|: 



Pn 



2r2 



Im{f,)H^ + Imik)Hj,)^o 



(39) 



and Hj^ have the same dependence on the functions as the radial sensitivity 

(my formulas are more compact but strictly 



functions and Hj^ 



m 



Tobie et al. 



2005b 



equivalent). For a satellite in synchronous rotation undergoing eccentricity tides, ^0 = 
21 eV5 (see Table [T| a nd n = w, in which case Pq is identical to the function htide{f) 
of Tobie et al. [2005b| (their Eq. (37)) except for a difference in sign explained below. 
After having shown that the global power is equal to An J htiiii.(r)r'^ dr , these authors 
correctly interpret htide{f) as the radial distribution of the dissipation rate per unit 
volume averaged over angles. 

Using variational principles, Tobie et al. [2005b prove that conservation of energy 
relates the sum of the strain and kinetic energies to the potential energy, the former two 
being integrated throughout the body whereas the latter is computable at the surface. 
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ECCENTRICITY TIDES: PATTERN A OBLIQUITY TIDES: PATTERN A 




-2.5 -2 -T5 -1 -0,5 0.5 1 ^ .5 2 2.5 



Figure 1: Basic patterns of tidal heating for eccentricity tides (left panels) and obliquity tides 
(right panels) for a body in synchronous rotation. In each case, the mean value has been sub- 
tracted from the angular function before normalizing it by its standard deviation. The 
origin of coordinates coincides with the sub-primary point when the primary is at pericentre 
(eccentricity tides) or at the ascending node (obliquity tides) . Patterns are repeated from 90° to 
-90°. 
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The dissipated part of this energy balance results in 



R 



Im{fi) + Im{K) ) dr 



5R 



Im{k2) 



(40) 



where R is the radius of the tidally perturbed body, k2 is the tidal gravity Love number 
of degree two and G is the gravitational constant. I obtain the correct sign in the right- 
hand side of Eq. (40) by correcting Eqs. (34)-(35) of Tobie et al. |2005b as follows: 
2/5 (-^s) must be replaced by its complex conjugate in their Eq. (34) so that the left-hand 
side of their Eq. (35) has an additional minus sign. 

The local power integrated over the volume of the body is thus given by 



E 



PdV 



An I P^r^dr 



5u;R 



Im{k2) {nRf • 



(41) 



This equation agrees with the formula that Zschau 1978 obtained for the energy dis- 



sipated during one tidal cycle by considering the dephasing between the tidal potential 
and the potential induced by tidal deformations (the formulation in terms of /m(/c2) was 
introduced by Platzman [1984] ). Therefore Eqs. (39), (40) and (41) prove the equivalence 
between the micro and macro approaches to tidal heating. 

Substituting the values of Table [l] into Eq. (41) and setting n = u, I get the well- 
known formula for the total dissipated power within a satellite in synchronous rotation 
\Cassen et al| |1980t [5^ate et al.\\l98^\Chyba et al.[|1989| , 



E 



^ , {ujRf /21 2 3.2, 



(42) 



which is valid up to second order in eccentricity e and obliquity / (terms proportional to 
esin/ vanish when averaged over angles). In the macro approach to tidal heating, the 
total power is often expressed in terms of the quantity 1/Q which Segatz et al. [1988 
define through 

Mfc2) = -^. (43) 

I inserted a minus sign in this equation because Q is usually assumed to be positive 
whereas Im{k2) is negative. There is however no simple relation between this l/Q and the 



specific dissipation function of the same name \Munk and MacDonald , 1960 , as discussed 



by Zschau [1978 and Segatz et al. 1988) . Contrary to what is sometimes said in the 



literature, Eq. (42) is valid when the body is inhomogeneous (but spherically stratified) 
and compressible. The macroscopic derivation of this equation does not even require a 
linear rheology. Note that existing extensions of Eq. ( |42[ ) to arbitrary eccentricity and 
obliquity | Wisdom , 2008 Levrard , 2008 assume that 1 /Q is proportional to frequency. 
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While this assumption makes computations easy, it does not correspond to any plausible 



rheology Efroimsky and Lainey , 2007 . Eq. ( 42 ) thus remains the best available formula 



for the total dissipated power within a satellite in synchronous rotation if the eccentricity 
and obliquity are not too large. 

In the Solar System, the obliquity of large satellites other than the Moon has not yet 
been observed except for Titan Stiles et al. , 2008 , [2010 but it is theoretically predicted 

Thus 



to be very small iGladman et al. 



Stiles et al. 


2008 


2010 


1996 


Peale 


1999 


Bills 



2005: Baland et a/. 2012 



obliquity tides are at present negligible for large satellites except the Moon where they 
contribute about 40% of the tidal heating. Heating patterns are however not observable 
on the Moon because radiogenic heating is dominant. Though Moon's obliquity may 
have been larger in the past, tidal heating due to obliquity tides was always smaller 
than radiogenic heating Peale and Cassen , 1978 Wisdom 2006 . Regarding exoplanets. 



Winn and Holman 2005 suggested that obliquity tides account for bloated 'hot Jupiters' 
but this idea has been criticized on the grounds that the corresponding Cassini state is 
unstable Fabrycky et aZ. , 2007 Levrard et al. . 2007 Peale, 2008. Tides could also 



be caused by forced librations in longitude Wisdom 



2004 



the effect of which can be 



computed by adding new terms proportional to P22 in the tidal potential (Eq. (88)). 



3.3 Surface flux 

Observations from space give indications on the distribution of surface heat flux but not 
on the dissipated power at depth. Computing angular variations of the surface heat flux 
thus involves assumptions about heat transport. Let me first define Xj ^ the fraction 
of the global power E due to the dissipation term J, 

R 

Im{fi) fj dr , (44) 

with a similar formula for Xk if J ~^ fx ~^ definition, J2j Xj + Xk ~ ^^ 

If the heat is radially transported to the surface, the surface heat flux due to tidal 
dissipation is given by 

T{e, <P) = {{XA + ^K)^ + XB^^+Xc^y (45) 
where J-q is the spatially-averaged heat flux at the surface, 

Under the assumption of radial heat transport, the coefficients Xj ^'^^ the weights for the 
three basic angular patterns of the surface flux. Non-radial transport mechanisms within 
the planet tend to average lateral variations in heating. On the other hand, an analysis 
including a more realistic modeling of heat transport within the planet (for example 
by convection) should also take into account lateral variations in viscosity due to tidal 
heating itself. I further discuss this topic in the Conclusions. Finally the heat passing 
through the crust can be channelled into heat pipes (for lo) or fractures (Enceladus' 
tiger stripes), thus further modifying the distribution of surface heat flux. 
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4 Interior structure and spatial patterns: toy models 



4.1 Thin shell 

4.1.1 Rheology constant with depth 

If the body contains an ocean beneath a thin icy shell (or an asthenosphere beneath 
a lithosphere), displacements do not vary much with depth within the shell and can 
be approximated by their value at the surface which depend on the displacement Love 
numbers h2 and l2- In the membrane limit of thin shell theory, these Love numbers are 
related by 

(47) 

where u is the viscoelastic counterpart of Poisson's ratio. This relation was previously 



derived for an elastic thin shell [Beuthe 2010, Eq. (D.12)]. Using the correspondence 
principle, I apply it here to a viscoelastic thin shell. The membrane limit of thin shell 
theory is appropriate for loads of large wavelength with respect to the shell th ickness d: 
the harmonic degree e of the load should satisfy £ <C l.S^/R/d \Beuthe\ |2008| Eq. (83)], 
that is d/R <^ 0.8 for tidal loads of degree two. This condition is satisfied in thin shell 
theory, in which one usually assumes that the shell thickness is less than 10% of the 
radius. 



When does the thin shell assumption break down for Eq. (47)? Consider an incom- 
pressible = I' = 1/2) two-layer body of uniform density in the liquid core limit for 
which the ratio /2/^2 is given by Eq. (114) with = 0. Deviations from the thin shell 
limit (^2/^2 = 3/11) are smaller than 10% if x = 1 — d/R > 0.9. This constraint coincides 
with the usual requirement of thin shell theory mentioned above and is satisfied in most 
models of icy satellites with a subsurface ocean. 



Eq. (47) is not very sensitive to rheology. For materials constituting planetary crusts. 



v ranges between 0.25 (silicates) and 0.325 (elastic ice. Gammon et al. [1983| ). The value 
of u depends on the rheological model. If the rheology is Maxwell, the real part of iy 
varies between the elastic value z/ (high viscosity) and the incompressible value 0.5 (low 
viscosity), while its imaginary part is always much smaller than its real part. Neglecting 
the imaginary part, we see that {1 + u) / (5 + changes from 1/4 to 3/11 (less than 10%) 
as P changes from 1/3 to 1/2. 



I compute weight functions in the thin shell approximation with Eqs. (25), (29) and 



(47): 



(/a ^ Jb ^ fc ^ Hfi) 



H 



K 



g- 

6 

11 



5 + 9 
- 29 



16 



1,0, 



9 11 

2 ' Y 



1 + 9 



H., 



(48) 



In absence of bulk dissipation, the relative weights of the patterns do not depend on 
the rheology: Pattern C dominates in the icy shell (fc/ ~ 3-8 for eccentricity tides) 
whatever the value of 9. Of course, the amount of tidal heating depends on the rheology 
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ECCENTRICITY TIDES: THIN SHELL 



OBLIQUITY TIDES: THIN SHELL 



90° 
60° 



-30° 



-90° 




-90" -60° -30" 0" 30" 60" 90° 
Longitude 




0° -60° -30° 0° 30° 60° 90° 
Lorgitude 



-2.5 -2 -1.5 -1 -0.5 0.5 1 1.5 2 2.5 



Figure 2: Pattern of tidal heating (^f^ + f^'c) due to eccentricity tides (left panel) or obliquity 
tides (right panel) within a thin compressible shell underlain by a liquid layer for a satellite in 
synchronous rotation (other details as in Fig. [I]). 



of ice through the common multiplying factor depending on /12 and D. Therefore, the 
dissipated power in a thin shell above a fluid layer (e.g. an icy crust above an ocean) is 
well approximated by 



8 io{nR)^ ^ 

^=3^:^ Mm) g2 



5 + z> 

where k is related to the ratio of bulk dissipation to shear dissipation: 

_ 11 Im{k)HK 



(49) 



(50) 



Fig. [2] shows the resulting spatial pattern for eccentricity tides and obliquity tides if 
there is no bulk dissipation (see Section 5.4 for a comparison with the literature). 



4.1.2 Love numbers 

In the formula for the dissipated power (Eq. (|49|)), the Love number h2 does not affect 
the spatial pattern but regulates the total amount of tidal heating. It can be computed 
with an interior model, for example the thin shell limit of the incompressible body of 



uniform density (see Eqs. (119)). I give here slightly more general formulas (obtained 



with the same method) valid for an incompressible body of nonuniform density having a 
fully rigid mantle surrounded by an ocean and a thin icy shell. The density and elastic 
structure below the mantle do not matter: a solid or liquid core of higher density may 
be present. The ocean and the icy shell have the same density. In this thin shell/rigid 
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mantle approximation, the Love numbers are given by 



h2 

h 
k2 



,(0) 



1 + 4°^ 



24 



d 



11 pgR R 



11 
5p 



h2 



(51) 



The parameters d, p and p are the sheh thickness, the density of the ice (or ocean) and 
the average density of the body, respectively. The factor h^^ is the value of /i2 in the 
limit of vanishing shell thickness: 



h. 



(0) 



5p 



bp — 3p 



(52) 



Eqs. (51) expanded at first order in d/R agree with Eqs. (3)-(6) of Wahr et al. 2006|. It 



is however important not to expand Eqs. (51) if one wants to compute good estimates 
of the imaginary part of Love numbers. Though Wahr et al. [2006 's formulas take 



into account the density difference between ice and ocean, they are not well-suited to a 
viscoelastic interpretation because Im{fl) is not necessarily of the same order as Re{p). 



Using Eqs. (49) and (51), I can check that the integration over the volume of the shell 



of the dissipated power yields the global power (Eq. (|42[)) as computed with the macro 
approach (equivalently, one can verify that Eq. ( [40| ) is satisfied). 

4.1.3 Depth-dependent rheology 

The above formulas implicitly assumed that the icy shell has a homogeneous rheology 
characterized by {p,v). Actually, the assumption of depth-independent rheology is not 
realistic. The viscosity of ice indeed depends on the local temperature T of the ice and 
thus varies a lot between the top (T ~ 100 K for Europa) and bottom (T ~ 273 K) 
of the icy shell, where it is at its melting point. This variation has a huge effect on 
the shear modulus /i, as can be seen for example in a Maxwell model \Wahr et al. 



2009 . However, Poisson's ratio P varies much less as I already argued in Section 



4.1.1 



In practice, one divides the shell into an outer conductive layer and an inner convective 



layer, the former being mainly elastic and the latter being viscoelastic iMcKinnon , 1999 



Dissipation mainly occurs in the convective layer whereas the conductive layer influences 
deformations. 

What is the effect of depth-dependent rheology on the thin shell formula for the 



dissipated power? There is no problem is letting Im{p) be depth-dependent in Eq. (49), 
but I2 is a constant so that (/i, v) cannot vary with depth in the formulas for Love 
numbers, Eqs. (47) and (51). In this last equation, we see that i) is constant [u = 1/2) 



and that jx is multiplied by the thickness of the shell d. In thin shell theory, it can be 



shown that the product pd in Eqs. (51 ) arises from the integration of the shear modulus 
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over the thickness of the shell \Beuthe 2008 : 



jld^ flavd 



jl dr . 



(53) 



shell 



In good approximation, Re{jlav) (resp. Im{jlav)) is determined by the conductive (resp. 
convective) layer where Re{jl) (resp. \Im{fl)\) is highest. Thus the real (resp. imaginary) 
part of the Love numbers, corresponding to deformation (resp. dissipation), is mainly 
determined by the conductive (resp. convective) layer. The substitution (1 — )• flav 



preserves energy conservation: one can either check that Eq. (40) is satisfied, or more 



explicitly that the integration over the volume of the shell of the dissipated power with 
depth-dependent fl yields the global dissipated power (Eq. ([42|)). Incompressibility (i) = 
1/2) is the simplest assumption for i) compatible both with depth-dependent rheology 



and with the absence of bulk dissipation. In contrast with Eq. (53), it would not make 



sense to simply replace u by its average in Eq. (47). Note that the averaging procedure 



described by Eq. (|53|) was heuristically proposed for an elastic icy shell by Wahr et al. 
[2006] 



In Section [5| I will compare thin shell predictions with the results of interior models 
of Europa and Titan. 



4.2 Homogeneous body 

The incompressible homogeneous body is the simplest interior model and serves as a 
useful approximation either (1) for bodies without liquid layers or (2) for a solid core 
surrounded by a weak layer (ocean or ice). The body is characterized by its radius R, 
density p and shear modulus /i. The problem of computing tidal deformations can be 
formulated in terms of the reduced radius f = r/R, the surface gravity g = {Att /'i)GpR 
and the reduced shear modulus defined by 



pgR 



(54) 



The solution to this venerable problem (going back to Lord Kelvin in the 19th century, 
see Love [1911| ) is for example given by Eqs. (7)-(9) of Peale and Cassen [1978| (in which 
{ae/r)'^V2mpq IS equivalent to in my notation). It is also derived in Appendix F (see 



Eq. (|107D). 

become 



After substitution of this solution, the weight functions given by Eqs. (24) 



(/a) /b) fc) 



\h2r 192 



g^ 25 



2(1 



' ) 



2 1 




(55) 



in which /12 is given by the first of Eqs. (105). While Hj^ 
the radial sensitivity function is given by 



(incompressible material), 



= — f2 (320 - 560 f2 + 259 f^) 
25 



(56) 
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Figure 3: Weight functions for the three basic spatial patterns of dissipation in a homogeneous 
incompressible body. The functions f j are the rescaled weight functions (Eqs. (30)), the sum 
of which is H^. The functions fg = 0.56 and — 0.84 are the equal- variance weight 

functions for eccentricity tides (see Eq. (38 1), which should be compared with = in order 
to determine which pattern dominates. 



with a maximum at r = 0.627. The comparison of the equal-variance weight functions 
fABC (^^^ ^^S- ^ ' shows that Pattern C dominates at all depths. 



Conservation of energy (Eq. (40)) can be analytically checked by radially integrating 
H and using the relations between /12, k2 and jl (see Eqs. (102) and (105)). I 



do 



not prove it here because it is a particular case of the two-layer model (see Eqs. (58) 
with z{0) = 19/2). Remarkably, weight functions depend on the rheology through the 
common factor |/i2p/(7^, so that the relative weights of heating patterns at a given radius 
do not depend on any physical parameter. The homogeneous solution is an example of 
the complete factorization of the overall deformation described by the tidal Love number 

h2. 

For an Earth-like elastic homogeneous body, allowing for compressibility changes /i2 
by 10 to 20% but k2 by only 2% [i.o^el|1911[ pp. 105-110]. If this result can be extended 
to a viscoelastic body, we expect sizable compressibility effects on the weight functions 
(depending on squared strains) but little change in the global power (depending on 
Im{k2))- Compressibility effects decrease with body size and thus do not affect much 



the radial sensitivity function for satellites of the solar system Tobie et al. 2005b 
Fig. 4]. The effect on weight functions remains however significant in small bodies with 



a liquid core (see Section 5.4). 



4.3 Incompressible two-layer model 

The radial functions m have analytical expressions depending on powers of r if the body 
is approximated as being (1) incompressible and (2) composed of spherically concentric 
homogeneous layers. Analytical formulas for the yi can be determined with the propaga- 



tor matrix technique [e.g. Sabadini and Vermeersen , 2004 and the help of mathematical 
software (I used Mathematica) though expressions are very lengthy except in the sim- 
plest cases. In this approach, one usually computes tidal deformations in the static limit 
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of infinite tidal period. 



4.3.1 Density effects 



Radial variations in density strongly affect the magnitude of the overall deformation. 
The relative weights of dissipation mechanisms, however, depend little on the density 
structure if the body is far from the fluid limit. In order to be more quantitative, let 
us consider a two-layer model with liquid core (radius Rc, density pc) and viscoelastic 



mantle (radius R, density pM, shear modulus /xm)- The rescaled weight functions (30) 
depend on three dimensionless parameters: the reduced core size x = Rc/R, the density 
contrast r/ = pc/Phi (supposed to be smaller than 5), and the reduced shear modulus 
of the mantle equal to p,red = P-M / {PqR) , with p the average density and g the surface 
gravity. I will assume that Im{p) <C Re{p), so that I can approximate p, and A by their 
real part when computing tidal displacements. I thus compute the weight functions in 
the elastic limit: the viscous response appears only in the formula for the power through 
the global factor Im{p). Density effects have a negligible effect on the rescaled weight 
functions if the body is far from the fluid limit, that is if 

l^redl >0.1(r7-l) . (57) 

If the density is uniform, the fluid limit is not relevant because the rescaled weight 
functions are independent of p^ed (see below). Fig.|4]shows how rescaled weight functions 
change as the density contrast between the liquid core {Rc = -R/2) and mantle increases 



from one to three. Density effects are small as long as \pred\ ^ 0.2 (see Eq. (57)), but are 



large if the body deforms like a fluid {pred — ^ 0). The condition (57) is generally satisfied 
in the dissipative layers of terrestrial planets and satellites of the solar system, though 
it could be violated for lo if the mantle is extremely soft with a shear modulus less than 



(ry — 1) GPa \Fischer and Spohn. 1990 Spohn, 1997|. It is true that the asthenosphere 



(if present) is characteristically much softer but its density contrast with the rest of the 
mantle is small. 



4.3.2 Uniform density and rheology 



Since the density structure has in most cases a small effect on the rescaled weights 
functions, I will from now on assume that the body is of uniform density. Rescaled 
weight functions in such models depend only on rheological contrasts between layers: 
the overall magnitude of the deformation factorizes (see Appendix F). In Appendix G, 
I give the analytical solution for tidal displacements within an incompressible body of 
uniform density p composed of a viscoelastic core (radius Rc, shear modulus pc) and a 
viscoelastic mantle (radius i?, shear modulus pm)- The model can be characterized by 
three dimensionless ratios (see Eqs. (108)): Rc/R-, Pc/fj-M and pred = Jj'M / [pgR)- The 



weight functions {fA^fB^fc) ^^"^ -^m '^^^ computed by substituting Eqs. (99), (109), 
(IT2|)-pT3|) into Eqs. ([^ and (Q . 

I will examine two interesting limits of the two-layer model, being (1) the dissipation 
within the mantle surrounding a liquid core, and (2) the dissipation within a soft layer 
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Figure 4: Influence of the density contrast rj — pc/pM between core and mantle on the rescaled 
weight functions ( |30[ ) for a two-layer incompressible body having a liquid core with Rc = 0.5R: 
rj — 1 (continuous curves), r] = 3 and pred — 0.2 (dashed curves), rj — 3 and pred — >■ (dotted 
curves). The densities pc and pM are adjusted so as to yield the same average density p and 
same surface gravity g; the relation between Pred and pm is thus not affected by the density 
contrast. 



above a rigid layer. These two special cases already display the different possible types 
of dissipation patterns. Besides they have the advantage that the rheology appears in 
the functions through the common factor /i2 describing the overall deformation. As 
a result, the relative weights of the three basic patterns depend only on the core radius. 
Dissipation patterns can thus be studied in these two models without specifying the 
rheology of the mantle. For an arbitrary core radius, the coefficients characterizing the 
functions in the mantle are given by Eqs. ( 109| ) and ( 112 )-( 113 ) in which the ratio ^ of 
shear moduli is either set to zero (liquid core) or tends to infinity (infinitely rigid core). 



4.3.3 Liquid core 

In the first limit case, the body consists of a viscoelastic mantle surrounding a liquid core 



of the same density. This model is equivalent to the one that Peale and Cassen 1978 



used to compute tidal dissipation in the Moon but my equations are simpler. If the core 



radius vanishes, the model is equivalent to the homogeneous case (see Eq. (107)). If the 



core radius tends to the surface radius, the model is a special case of the thin shell model 



discussed in Section 4.1 (compare Eqs. (119) with Eqs. (47) and (51)). 

Figs.[5]^a,b,c,d) shows the rescaled weight functions in the mantle surrounding a liquid 
core for different values of the core radius. The location of maximum average dissipation 
(maximum H^) shifts from the value r/R = 0.627 (homogeneous model) to the core- 
mantle boundary as the core size increases. The comparison of the weight functions shows 
that Pattern C is generally dominant at all depths. At the surface, Pattern A slightly 
dominates Pattern C if Rc ~ 0.5 — 0.75R (Figs. [5](b,c)), in which case the resulting 
pattern is mostly of degree four (the component of harmonic degree two approximately 



cancels, see Eqs. (36)) but with a pattern reversed with respect to Pattern B, as in 
Fig. 6(f) of Tobie et a/.| [2005b] . Pattern A is however negligible at greater depth so that 
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the radially-integrated power is dominated in all cases by Pattern C. The contribution 
of Pattern B is significant if the core is small (/^//(^ maximum for Rq ~ The 
influence of core size is summarized in Fig. ^a) showing the relative contribution of each 
dissipation term to the total power. 



4.3.4 Rigid core 



In the second limit case, the body consists of a viscoelastic mantle surrounding an 
infinitely rigid core of the same density, though 'mantle' and 'core' are only names for 
the layers and do not necessarily refer to their physical counterparts. This toy model is 
relevant to at least two realistic configurations. First, the model is a good approximation 
for dissipation within an icy layer just above a rocky mantle, as ice is much softer than 
rocks. Second, the model describes the essentials of asthenospheric tidal heating in which 
tidal dissipation mainly occurs within a soft layer at the top of the mantle. 

Figs. [5][^e,f,g,h) show the rescaled weight functions in the mantle surrounding an 
infinitely rigid core for different values of the core radius. The comparison of the weight 
functions shows that Pattern B dominates near the core-mantle boundary where shear 
is maximum and decreases to zero near the surface. 

If the soft layer is thick, the weights of the different patterns are of comparable 
magnitude, with Pattern C, A or B successively dominating as the thickness of the 
rigid layer increases from zero to the surface radius (see Fig. [6][^b)). For example, a 
soft layer with thickness R/2 leads to weights satisfying on average fs ^ fc ^ ^° 



that Pattern A dominates by virtue of Eq. (37). As the soft layer becomes thinner (say 
^ 0.75i?), Pattern B dominates Pattern C in the lower part of the layer with the 
reverse situation near the surface. Therefore, the dissipation pattern within an icy layer 
in contact with the mantle is given by Pattern B at the bottom of the layer and by a 



mix of Patterns A and C at the top of the layer, as in Figs. 10(a,b,e,f) of Tobie et al 
|2005b| . 

When the soft layer is very thin (say Rc ^ 0.9-R), Patterns A and C become negligible 
because dissipation is mainly caused by shear heating due to the coupling of the soft 
layer with the rigid layer. This mechanism accounts for the dominance of Pattern B in 
models of asthenospheric dissipation in lo \Segatz et al. 1988 . Such models include a 
rigid lithosphere above the asthenosphere so that shear dissipation also occurs at the 
asthenosphere-lithosphere boundary, yielding a U-shape for the curve /b (see application 



the sub- Jovian and anti- Jovian points \Lopes-Gautier et al., 


1999 


Tackley et al. 


2001 


Kirchoffet al., 2011 


Veeder et al. 


2012 


Hamilton et al. 


2012 


, models of asthenospheric 



dissipation (less heating at the poles) have been more popular than models of mantle 
dissipation (more heating at the poles) . There is however no evidence for the degree- four 
structure exhibited by Pattern B. Another possibility is provided by a model having a 
rigid lower mantle and a soft upper mantle of large thickness, in which case Pattern A 
gives a significant contribution partly cancelling the degree-four term due to Pattern B 
and leading to a degree-two structure reverse of Pattern C (for example + "^^b ~ 
5^Q — from Eq. (37)). The influence of asthenosphere thickness is summarized in 
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Figure 5: Rescaled weight functions (30) for the three basic spatial patterns of dissipation in the 
mantle of a two-layer incompressible body with either a liquid core (left panels) or an infinitely 
rigid core (right panels) . Each row corresponds to a different core radius with Re /R taking the 
values 0.25, 0.5, 0.75, 0.9 from top to bottom. The big dots on the right-hand side of Panel (d) 
are the values of the functions in the thin shell limit {Rc ~^ R). Note the much larger scale 
in Panels (g) and (h). Equal- variance weight functions are not shown for clarity; they could be 
drawn as in Fig. [s] from the relations valid for eccentricity tides: = 0.56/^ and — 0.84/, 



c- 
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(a) Liquid core (b) Rigid core 
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Figure 6: Contribution of Pattern 



layer incompressible body of uniform density as a function of the core radius 
infinitely rigid core. 



(see Eq. (44|) to the global dissipated power for a two- 

(a) liquid core, (b) 



Fig. [6Fb) showing the relative contribution of each dissipation term to the total power. 



4.3.5 Conservation of energy 

In the two limit cases (liquid or rigid core of reduced radius x) , dissipation occurs only in 



the mantle. It is thus not difficult to analytically check conservation of energy (Eq. (40)) 
by noting that 



Im{jl) 



|/^2(X)|- 



5 pgR Im{h2{x)) 

2 4x) |/l2(x)|2 



(58) 



where h2{x) and z{x) are defined by Eqs. (116)-(118) and /12 = (5/3)A;2 (see Eq. ( 102[ )). 
Considering still the two limit cases, I can use either one of Eqs. (58) to write a compact 



formula for the global power dissipated by eccentricity tides in terms of core size and 
mantle rheology: 



E 



63 2 i^^RY 



Im.{fi. 



red) 



z{x) 



(59) 



4 G " |-|^ _^ z{x)flred\ 

This equation serves to estimate the core size and the rheology yielding the required 



global power output (e.g. Fig. 1 in Peale et al. 1979 and Fig. 1 in Cassen et al. 1980| ). 



5 Real bodies: Europa, Titan and lo 

5.1 General features of interior models 

I will show here that the toy models of Section |4] predict well the spatial patterns of 
tidal heating within each layer of a real body that can be approximated by concentric 
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Table 4: Global parameters for Europa, Titan and lo: surface radius i?, average density p, 
normalized axial moment of inertia Mol=C/MR^, surface gravity g and tidal angular frequency 
ui. 





R (km) 


p (kg/m3) 


Mol 


g (m/s2) 


uj (rad/s) 


Europa 


1561 


3013 


0.346 


1.31 


2.05 X 10-5 


Titan 


2575 


1881 


0.341 


1.35 


4.56 X 10^*^ 


lo 


1822 


3527 


0.378 


1.80 


4.11 X 10^5 



homogeneous layers. Toy models however do not predict the magnitude of dissipation in 
the layer, this factor depending on the complete structure of the body. Thus comparisons 
cannot be done between layers unless the tidal response of the real body has been 
computed. 

Choosing as illustrations the satellites Europa, Titan and lo, I approximate their 
internal structure with four or five layers that are homogeneous and incompressible. 
In its general features, the interior model for Titan is also relevant to Ganymede and 
Callisto, which probably have a subsurface ocean sandwiched between thick ice layers 

Triton is yet another candidate 



Spohn and Schubert, 2003 



et al. 



Hussmann et al. 



2006 



Gaeman et al. , 2012 . Global physical parameters (see Table [4|) are taken from Schubert 



2009 and Sotin et al. 2009 for Europa, from less et al. 2010| for Titan and 



from Moore et al. 2007| for lo. Solid layers are viscoelastic with Maxwell rheology 
(see Eq. (74)). Each layer is characterized by its thickness 6R, density p, elastic shear 
modulus /i, and viscosity rj. As in Section 4.3, I compute the functions with the 



propagator matrix technique for an incompressible body in the static limit. 



5.2 Interior models of Europa, Titan and lo 

I model Europa with five layers: a liquid iron core, a silicate mantle, an ocean and an 
icy crust stratified into convective and conductive layers. Table [5] gives the values of 
the internal physical parameters; they are based on those of Tobie et al. 2003 Table 1 



and Fig. 5], except that I adjusted the core and mantle radii so as to obtain the correct 
mass and moment of inertia. Constraints on density structures and on rheology are 
reviewed by Schubert et al. [2009 and Sotin et al. [2009 , respectively. Estimates of the 
ice thickness range from 100 m for the upper elastic layer to 30-50 km for the total ice 
thickness (see reviews by Billings and Kattenhorn [2005 and Nimmo and Manga 2009| ). 
I model Titan with five layers: a rocky core, a mantle of high-pressure ice, an ocean 



and an ice layer stratified into convective and conductive layers. Fortes 2012 proposed 



several internal structures, among which I choose his 'dense-ocean' model (see Table [g]), 
which yields the correct moment of inertia and is compatible with the value of tidal k2 
measured with Cassini less et al., 2010 2012 . For the rocky core, I assume the same 



rheology as in Europa's mantle. The elastic shear modulus of ice is determined from 
pVg, with the S-wave velocity given by = 1880 ms" 



Tobie et al. 



2005a 



while the 

ice viscosity is chosen to be the same as in Europa's crust. It is not clear whether the 
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Tabic 5: Interior model of Europa 



Layer 


SR (km) 


P (kg/m^) 


H (GPa) 


rj (Pa.s) 


Iron core 


587 


5500 








Silicate mantle 


840 


3500 


70 


1020 


Ocean 


114 


1000 








Convective ice 


12 


920 


3.3 


10" 


Conductive ice 


8 


920 


3.3 


1020 



Table 6: Interior model of Titan 



Layer 


5R (km) 


P (kg/m^) 


(GPa) 


rj (Pa.s) 


Rocky core 


1984 


2650 


70 


1020 


Icy mantle 


241 


1351 


4.8 


10" 


Ocean 


130 


1000 








Convective ice 


10 


931 


3.3 


10" 


Conductive ice 


90 


931 


3.3 


1020 



icy crust convects [e.g. 



Sohl et al. 


2003 


Mitri and Showman. 


2008 


Nimmo and Bills 



2010 



I thus assume that convection is limited to the bottom 10% of the icy shell in 
order to illustrate the effect of a thick conductive ice layer above an ocean. This interior 
model yields Re{k2) = 0.57. 

I model lo with four layers: a liquid iron core, a silicate mantle of uniform density 
but rheologically divided into a deep mantle and an upper asthenosphere, and a crust 
(or lithosphere). Table [7] gives the values of the internal physical parameters. Given 
plausible values for the crust thickness and density, I constrain the density structure 
with the average mass and moment of inertia \Moore et al. 2007| . As there is still 
one free parameter, I set to 0.41 the ratio of core radius to deep mantle radius, the 
same value as in Europa's interior model. With this choice, I intend to illustrate the 
predictive power of the toy models of Section [4| apart from their overall magnitude, 
weight functions will look the same in Europa's mantle and lo's deep mantle. 

I set the asthenosphere thickness to 50 km which is the lower bound for the conduc- 
tive layer detected with Galileo \Khurana et al. 2011 . Since lo's internal temperature 



distribution and degree of partial melting are unknown \Keszthelyi et al. , [2007 \Moore 



et al, 2007; Kohlstedt and Mackwell . 2009|, the viscoelastic parameters {iJ,,r]) are not 



determined from known properties of rocks but are chosen so as to satisfy the global 
constraint of Im{k2) = — |^2|/Q = —0.015 it 0.003 Lainey et al. . 2009] . This con- 
straint comes from astrometric observations of the secular accelerations of lo, Europa 
and Ganymede. The substitution of Im{k2) in Eq. (42) with e 
a total power of = (9ib2) x 10^'^ W compatible with the estimates 6 
heat flux measurements 



0.0041 and / = yields 
16 X 10^3 w fj-om 



Moore et al. , 2007 . For the deep mantle, I adopt (/i, rj) values 



suggested by Spohn 1997|. For the asthenosphere, I choose (/i,?/) so as to satisfy the 
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Table 7: Interior model of lo 



Layer 


SR (km) 


p (kg/m^) 


fi (GPa) 


r/ (Pa.s) 


Iron core 


716 


6767 








Deep mantle 


1026 


3349 


10 


IQis 


Asthenosphere 


50 


3349 


5 X IQ-^ 


2.5 X 10^° 


Crust 


30 


2750 


70 


1020 



global constraint on dissipation. In this model, tidal heating is nearly equipartitioned 
between deep mantle (49%) and asthenosphere (51%), the dissipation in the lithosphere 
being always negligible. 

5.3 Dissipation weight functions 

The spatial patterns resulting from weight functions have been discussed at length in 
Section |4j I thus restrict myself here to comparing the weight functions within Europa, 
Titan and lo with the weight functions of the toy models of Section |4j 

Let me start by examining dissipation weight functions in Europa's mantle, Titan's 
rocky core and lo's deep mantle. In Europa's mantle (Fig. [7][^a)), weight functions look 
like those in the mantle of a two-layer body with liquid core (the curves are intermediate 
between Figs. [5]^a) and[5]^b) as the core radius is 41% of the mantle radius). This sim- 
ilarity could be expected, because the free-slip condition at the mantle/ocean interface 
mimics one of the boundary conditions at the surface of the body. The presence of an 
ocean has also the effect of greatly reducing the amplitude of weight functions in the 
mantle because deformations of the mantle are due to gravitational coupling: they van- 
ish if the core, mantle and ocean have the same density. The surrounding ocean actually 
provides most of the tidal response of the satellite. In Titan's core (Fig. [7][^b)), weight 
functions look like those in a homogeneous body (compare to Fig. [s]), which is not sur- 
prising since the icy mantle does not impose a strong shear coupling at the core/mantle 
boundary. In lo's deep mantle (Fig. [Tjjc)), weight functions look like those in Europa's 
mantle though their magnitude is one order of magnitude larger. This resemblance is 
explained by similar boundary conditions: on the one hand the mantle encloses a liquid 
core with the same ratio of core radius to mantle radius (41%), and on the other the 
mantle is surrounded by a much weaker layer (water for Europa, asthenosphere for lo). 

Next, I compare dissipation weight functions in the crusts of Europa, Titan and lo. 
In Europa's icy crust (Fig. [Tj^d)), weight functions are nearly constant as is expected in 
a thin shell (compare to Fig. [sj^d) where the shell is much thicker). In particular, the 
stratification of the ice layer has nearly no effect on the weight functions (a zoom on 
reveals a discontinuity at the convective/conductive transition, but this weight function 



is close to zero in the icy shell). The thin shell/rigid mantle approximation (Eqs. (48) 



and ( 51 )-( 53 )) predicts well the values of the weight functions in the icy shell (big dots in 
Fig.[7p)). One measure of comparison is Im{k2), as it summarizes the effect of internal 



structure on the global dissipated power (see Eq. (42)). The interior model of Europa 
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specified in Table [s] yields /m(/c2) = —3.58 x 10~'^ whereas the thin shell/rigid mantle 
approximation gives Im{k2) = —3.53 x 10~^, a difference of less than 2% which is mainly 
due to the elasticity of the mantle and the fluidity of the core (the dissipation in the 
mantle is negligible). 

In Titan's icy crust (Fig. [Tj^e)), the dominant weight function decreases by 20% 
between the bottom and top of the shell; this variation is compensated by an increase 
in so that is approximately constant with depth. Thin shell formulas thus do 
not yield as good results as for Europa, but are still approximately valid given the large 
uncertainties in the physical parameters. The thin shell/rigid mantle approximation for 
the weight functions (big dots in Fig. [7](e)) is off by about 10% mainly because of the 
large density difference between ocean and icy shell in the 'dense ocean' model. In the 
thin shell approximation with /12 given by the interior model (triangles in Fig. [7|^e)), 
estimates of weight functions are off by only a few percents from their average values 
over the shell thickness. Spatial patterns in Titan's icy shell are thus well predicted by 
this type of thin shell approximation. As for Europa, the stratification of the ice layer 
has nearly no effect on the weight functions. Regarding global dissipation, the interior 
model yields Im{k2) = —7.9 x 10~^ whereas the thin shell/rigid mantle approximation 
gives Im{k2) = —8.3 x 10^'^, a difference in the imaginary part of less than 5% (the 
overestimated weight functions are partially compensated by the neglected dissipation 
in the icy mantle, which amounts to 4% of the total). 

In lo's crust, weight functions are nearly constant and are well predicted by the 
thin shell approximation with /i2 given by the interior model (triangles in Fig. [7][|f)). 
The weakness of the asthenosphere underlying the crust results in a nearly free-slip 
condition for the crust and accounts for the success of this approximation. The elasticity 
of the asthenosphere cannot however be neglected when computing the magnitude of the 
deformation. If the asthenosphere is absent, the crust is strongly coupled to the mantle 
and the thin shell approximation cannot be used. In that case, the values of the weight 
functions in the crust are well approximated by their values at the top of the mantle, 
which are sensitive to the size of the liquid core (see Figs. [5][^a,b,c,d)). For example, 
/a ~ fc ratio of core radius to mantle radius is about one-half (/^ is always close 

to zero in the crust). 

Finally, I compare dissipation weight functions in Titan's icy mantle and lo's astheno- 
sphere. In Titan's icy mantle (Fig. [sj^a)), weight functions resemble those in the mantle 
of a two-layer body with an infinitely rigid core of radius Rc = 0.9i2 (see Fig. ^h)). 
This result is explained by similar boundary conditions: free-slip condition at the man- 
tle/ocean interface and rocky core much more rigid than the icy mantle. Besides, the 
thickness of Titan's icy mantle is about 10% of the mantle radius as in Fig.[5]^h). In the 
lower half of lo's asthenosphere (Fig. [sj^b)), weight functions look like those in Titan's 
icy mantle (or equivalently to those in a thin weak layer above a rigid core) but with 
an even stronger dominance of The thickness of lo's asthenosphere in the chosen 
model is indeed less than 3% of the total radius, which is much less than in Fig. [5](h). In 
the upper half of lo's asthenosphere, weight functions are nearly the mirror image of the 
functions in the lower half because the upper boundary is not free but instead in strong 
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shear with the much more rigid crust. 



5.4 Comparison with the literature 



The formahsm developed in this paper is consistent with previous results found in the 



literature. First, the formula for the global power (Eqs. (41)-(42)) obtained by integrat- 



ing the local power over the volume of the body agrees with the formulas that Zschau 
[1978] obtained using the macro approach to tidal heating. Second, the formula for the 



radial distribution of the dissipation rate (Eq. (39)) is the same as the one derived by 

are the same as 



Tobie et al. 2005b 



et al 



Tobie 



In particular, the functions and Hj^ 
2005b| s radial sensitivity functions with the same name. These equivalences have 



already been discussed in Section |3.2[ Third, the dominance of Pattern C obtained in 
Section [42] in a homogeneous body agrees with the results of Peale and Cassen [1978 
Figs. 1-2] both for eccentricity tides and for obliquity tides. 

Consider now the two maps of lo's surface flux that Segatz et al. [1988| computed for 
mantle dissipation (Model A) and asthenospheric dissipation (Model B). These authors 
compute tidal strains in the incompressible limit with the propagation matrix method. 
It is difficult to reproduce exactly their figures because of some inconsistencies and 



missing information. Spohn 1997 and Schubert et al. 2004 already noted that Segatz 



et al. 1988 forgot a factor 2tt in the angular frequency w (the error is limited to the 
computation of the complex shear modulus and neither affects the factor u multiplying 
Im{fl) nor from the tidal potential). Thus all viscosity values in Segatz et al. 



1988 



should be divided by 2-it since the Maxwell viscosity is always multiplied by oj. Besides 
these authors do not not clearly state the viscosity used for the mantle in Model A and 
for the asthenosphere in Model B. As regards Model A, I duplicate the map of surface 

if the rheology of the mantle is specified by 



1988 



flux shown in Fig. 8 of Segatz et al. 
/i = 10^*^ Pa and r] = 3/(27r) x 10^^ Pa.s (as in their 'preferred' model shown in their Fig. 3, 
though the viscosity mentioned in their text is two-thirds of this value) . My choice leads 
to a sub-Jovian flux of 0.6 Wm~^ as stated in their paper; the corresponding total power 
is 10% smaller than their constraint of = 6 x 10^^ W, with |/c2| = 0.246 and Q = 28 
(note that the values \k2\ = 0.25 and Q = 36 mentioned by Segatz et al. [1988 yield 
a power 28% lower than their power constraint). The patterns (^^,^'^,^(^) have the 
weights (7.2, 22.3, 70.5)% in the surface flux. As regards Model B, I duplicate the map of 



surface flux shown in Fig. 10 of Segatz et al. [1988 if the rheology of the asthenosphere 
is specified by ;U = 10^ Pa and rj = 1.5/(27r) xl(? Pa.s, in which case the total power 
constraint is satisfied (|/c2| = 0.82 and Q = 84) and the maximum flux reaches 2.4 Wm~^ 
as stated in their paper. If I use the viscosity that Segatz et al. 1988] claim to have 
used for their Fig. 5 {t] = l/(27r) x 10® Pa.s), I obtain a power 50% lower than their total 
power constraint. Whatever the precise value of the viscosity, the distribution of surface 
heat flux is nearly 100% Pattern B. 



What is the correspondence between the thin shell formulas of Section 4.1 and the 

? 



results of\Ojakangas and Stevenson\ [1989a 



The substitution of Eq. (47) into Eqs. (26) 



reproduces in a simpler way and with more generality (i.e. allowing for compress- 
ibility) the strains computed by these authors (their Eqs. (B.20)-(B.25)), with their first 
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Figure 7: Weight functions in the mantle/rocky core/deep mantle (left panels) and in the crust 
(right panels) of Europa, Titan and lo. Interior models are specified in Tables [5] [6] and [Tj Note 
the change of scale on y-axis between panels. In the right panels, the vertical line separates the 
convective layer (if present) from the conductive layer. The big dots are the values of weight 
functions in the thin shell/rigid mantle approximation used for floating icy crusts (|ft.2|=l-21/1.54 
for Europa/Titan), whereas the triangles are the values of the weight functions in the thin shell 
approximation with \h2\ given by the interior model (|/i2|=l-21/1.45/0.99 for Europa/Titan/Io). 
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(a) Titan: high-pressure ice mantle (b) lo: asthenosphere 




radius [km] radius [km] 



Figure 8: Weight functions in (a) the high-pressure ice mantle of Titan and (b) in the astheno- 
sphere of lo. Interior models are specified in Tables [6] and [7] The scale is not the same in the 
two panels. Curves and Jq cannot be distinguished in Panel (b) because of the large scale. 



three coefficients gi (the coefficient 54 is not needed) corresponding to 



(91,92,93) = e{ry[,y^,y^] 



(60) 



In the right-hand side, one must set P = 1/2 (incompressible limit) and /i2 = (5/2)C/j 



[Cr = 0.5 corresponds to the choice p = 3p in Eq. (52), this density ratio being appro- 



priate for Europa). Besides, I can duplicate Fig. 1 of Ojakangas and Stevenson 1989a 



showing the quantity ef- in units of (uj-yeCu) /2, in which e is the eccentricity and 



7 = CO R/g. Combining Eq. (49) and Eq. (73), I obtain the correspondence 



-2 



9 



c 



(61) 



Lastly, consider the spatial patterns of dissipation in icy satellites computed hy \Tobie 
et al. [2005b who took into account compressibility effects. First, these authors examine 
the dissipation within the rocky mantle which is either homogeneous or includes a liquid 
core, the liquid core being either small or large. In their Fig. 6, they show the spatial 
patterns of maximum dissipation rate, surface dissipation rate and surface heat flux 
(if radial heat transport) for these three models. As regards the homogeneous model, 
I can duplicate Figs. 6(a,d,g) of Tobie et al. 2005b using my Eqs. (55) valid in the 



incompressible limit, because compressible effects are very small for a homogeneous body 
comparable in size to icy satellites. As regards the two models with liquid core, | Tobie\ 



et al. [2005b] already observed in their Fig. 4 that finite compressibility increases (resp. 
decreases) dissipation at the bottom (resp. top) of the mantle, but that the two effects 
approximately cancel when integrating on the radius. In particular, compressibility leads 
to the damping of the weight function at the lower and upper mantle boundaries and 
to an enhancement of f^-, in the lower part of the mantle in comparison with my Fig. [5|b) 
(note that the effect of compressibility on weight functions can vary a lot depending on 
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the value of the reduced shear modulus and the type of boundary conditions). I can 
duplicate Figs. 6(b,e,h) and Figs. 6(c,f,i) of Tobie et al. 2005b with the parameters 



given in their paper by computing the functions with a Fortran code allowing for 
compressibility which was originally developed by Dahlen for the Earth Dahlen, 1976 



(another version of this code is used by Wahr et al. [2006 , 2009] ) . The basic differences 
between the resulting patterns are well explained by the toy model of a viscoelastic 
mantle above a liquid core. For example, the component of harmonic degree four is 
much more visible at the top of the mantle when the core is large {Rc ~ -R/2, Fig. 6(f) 
of Tobie et al. [2005b| ) than when the core is small {Rc 
|2005b| 



because in the former case f. 



Section 4.3.3). 



R/3, Fig. 6(e) of \Tobie et al.\ 
Jq at the top of the mantle (see discussion in 



Next ^obie et a/.| [2005b examine the dissipation within icy layers for three different 
models: a thick icy layer in contact with the rocky mantle, a thin icy layer above a sub- 
surface ocean and a thick high-pressure ice layer sandwiched between the rocky mantle 
and a subsurface ocean. In their Fig. 10, they show the spatial patterns of dissipation 
at the top and bottom of these icy layers. I choose to reproduce their results with the 
incompressible propagator matrix method because this method works well when there 
are large rheological contrasts between layers. I can duplicate nearly all spatial patterns 
of Tobie et al. |2005b| with the parameters given in their paper: their Figs. 10(a,e) 



result from Pattern B only (bottom of strongly coupled icy shell), their Figs. 10(b,f) 
result from a mix of Patterns A and C (top of strongly coupled icy shell) whereas their 
Fig. 10(d) corresponds to the pattern in a thin icy shell above an ocean. These results are 
already well predicted by the toy model of a viscoelastic mantle above a rigid core (see 



discussion in Section 4.3.4) and are also similar to what I obtained for realistic models 



of Europa and Titan (see my Figs. [7|jd,e) andjsj^a)). The pattern in their Fig. 10(c) is 
not absolutely correct. Showing the pattern at the bottom of the thin icy shell in their 
'Europa-with-ocean' model, it should be the same as the pattern at the top of the shell 
(their Fig. 10(d)) because the ice shell is extremely thin (1 km is less than 0.1% of the 
radius). In that case, weight functions are nearly constant through the shell (see my 
Figs. [7]^d,e)). Compressibility does not play a role since its influence on spatial patterns 
is negligible in a thin shell. The difference between Figs. 10(c) and 10(d) of Tobie et al. 
[2005b| is however small. 



6 Conclusions 

In this paper, I showed that three basic spatial patterns are enough to describe the 
angular dependence of the power dissipated by tides in a spherically stratified body. 
This means that the dissipated power at a given radius is the sum of three terms, 
each term being the product of a radial function (or weight function) depending on the 
internal structure and an angular function (or spatial pattern) depending only on the 



tidal potential (Eqs. (22)-(24)). One term results from radial-tangential shear strains, 



another from tangential strains while the remaining term (usually the subdominant one) 
receives contributions from radial strains, tangential strains and volume change. This 
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decomposition presents the advantage that the 3D problem of predicting the spatial 
pattern of dissipation at all radii is reduced to the ID problem of computing the tidal 
displacements. As a further simplification, the three basic angular functions can be 
expressed as linear combinations of the harmonic components of degrees zero, two and 
four of the squared norm of the Fourier-transformed tidal potential (Eqs. (36)). The 
three basic patterns of tidal dissipation (Fig. [T]) are thus known once the tidal potential 
is specified, the most important case being a body in synchronous rotation undergoing 
eccentricity tides. 

The angular average of the spatially-dependent dissipated power yields the radial 
distribution of the dissipated power (already known from conservation of energy), the 
integration of which gives the classical formula for the global dissipated power. This 
last quantity is independently known, being proportional to the surface integral of the 
product of the tidal potential with the induced potential. The formula for the dissipated 
power derived in this paper thus completes the connection between the microscopic 
approach to tidal dissipation, based on the computation of strains, and the macroscopic 
approach, based on the phase delay of the induced potential (or on the phase delay of 
the tidal bulge). 

I studied the effect of internal structure on the dissipation pattern with different 
types of toy models. In a thin icy shell above an ocean, the dissipation pattern is 
nearly invariant through the shell thickness while the dissipation amplitude varies with 
depth as the rheology changes from elastic (conductive layer) to viscoelastic (convective 
layer). In an incompressible homogeneous body (also relevant to the solid core), the 
dissipation pattern depends on depth but not on rheology. The third toy model is 
the incompressible body with two viscoelastic layers. I demonstrated that the density 
contrast between the two layers weakly affects weight functions ratios. More generally, 
assuming uniform density in a multi-layer body leads to the factorization of the overall 
deformation (described by the Love number /12), so that the spatial pattern in each 
homogeneous layer depends only on contrasts of shear moduli at the layer boundaries. 

An important subset of the two-layer model consists of bodies with a liquid core. 
I showed how the weight functions change as the core radius increases from zero (ho- 
mogeneous model) to the surface (thin shell model). In all cases, the same pattern is 
dominant, resulting in more heating at the poles and less heating at the equator, at least 
if the body is in synchronous rotation and undergoes eccentricity tides (this result was 
already known). Adding a liquid layer and an icy shell above the viscoelastic mantle 
reduces the magnitude of weight functions within the mantle but does not change their 
ratios, so that the spatial pattern within the mantle remains approximately the same. 

Another subset of the two-layer model consists of bodies with a viscoelastic mantle 
surrounding a core of infinite rigidity. While this model may seem at first sight unphys- 
ical, it describes well the increase in radial-tangential shear dissipation occurring at the 
interface between a soft and a rigid layer. This situation could for example occur in 
icy satellites if the icy layer is in contact with the rocky mantle. This behavior is also 
typical of the asthenospheric dissipation that could explain lo's volcanism, except that 
radial-tangential shear dissipation is also large at the top of the asthenosphere where 
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it is in contact with the more rigid hthosphere. If the soft layer is thin, the spatial 
pattern for eccentricity tides within a body in synchronous rotation corresponds to less 
heating at the poles and more heating in the equatorial zone, with a significant com- 
ponent of harmonic degree four (the maxima are thus at the north and south of the 
sub- and anti-primary points). This pattern was already known from internal models of 
lo. If the soft layer is thick, the three basic patterns combine with comparable weights, 
resulting in a new type of dissipation pattern, with maximum heating at the equator as 
in asthenospheric dissipation but with a lower content of harmonic degree four. 

In summary, there are three classes of spatial patterns associated with 1) mantle 
dissipation (with the icy shell above an ocean as a particular case), 2) dissipation in a 
thin soft layer (asthenosphere or icy layer in contact with the mantle or core) and 3) 
dissipation in a thick soft layer (deep asthenosphere or thick icy layer above a solid core) . 
I demonstrated the generality of the inferences drawn from toy models by computing 
weight functions for realistic internal structures of Europa, Titan and lo. Finally, I 
duplicated spatial patterns previously obtained in the literature. 

The formalism described in this paper rests on two crucial assumptions: linear rhe- 
ology and spherical symmetry of the internal structure. Further assumptions on heat 
transport are necessary in order to compute the surface heat flux. The easy way out is 
to assume radial heat transport since it directly leads to the linear superposition of dis- 
sipation patterns at all depths whereas lateral heat transport tends to average patterns 



in complex ways. For example, Tackley et al. 20011 computed 3D convection within 



lo using as input tidal heating predicted by spherically stratified models and neglecting 
lateral variations of viscosity. At low Rayleigh number, the heat flux follows the input 
dissipation pattern but small-scale instabilities spread the heat flux as the convective 
vigor increases. In these simulations, the Rayleigh number is still much lower than pos- 
sible values within lo. In another paper, Tackley 2001 used 2D simulations at a more 
realistic Rayleigh number in order to show that large-scale lateral flows tend to average 
the dissipation pattern. The averaging effect due to convective heat transport is however 
counteracted by the feedback of nonuniform tidal heating on local viscosity. Using codes 
coupling convection and tidal heating, Behounkovd et al] [2010] and Han and Showman 



[2010| showed that tidal heating is significantly higher (resp. lower) in warm (resp 



colder) convective plumes. With their 3D code, Behounkovd et al. |2010| obtained dis- 



sipation patterns in Enceladus' mantle roughly similar to the pattern predicted by the 



method assuming internal spherical symmetry but with more small-scale structure ( Han 



and Showman [2010 's code is 2D and cannot predict global patterns). Computations 



coupling convection and tidal heating in 3D are very complex and have neither been 
applied to all types of dissipation (such as asthenospheric heating) nor been system- 
atically compared to the results obtained with the method assuming internal spherical 
symmetry. For these reasons I believe that the analytical approach presented here will 
remain a useful and practical tool to predict tidal dissipation patterns. 
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Appendix A. Dissipated power from strains 

I show here how the dissipated power can be expressed in terms of viscoelastic strains 
in the Fourier domain. The formulas are weh-known but not always correctly stated in 
the literature. Tidal deformations induce stresses aij{t,x) and strains €ij{t,x.) at each 
point X within the planet. If tides operate at only one angular frequency uj, stresses and 
strains can be written in terms of complex Fourier components cjjj and iij: 



aij{t, x) = - aij exp{iujt) + c.c. , 
eij {t, x) = ^ iij exp(ia;t) + c.c. , 



(62) 



where c.c. means complex conjugate. The dependence on x of the Fourier components is 
implicit. If there are several tidal frequencies, the Fourier decomposition involves a sum 
over frequencies. The quantities aij and eij (or dij and iij) are the physical components 
of the stress and strain tensors, that is the Cartesian components in a local system of 
rectangular Cartesian coordinates with axes tangent to the coordinate curves at the 
point \Fung 1965 Malvern, 1969 . In spherical coordinates, this means that the basis 
vectors are normalized. 

In the Fourier domain, the constitutive equations for an elastic solid take the form 
\Fung 



1965 Malvern, 1969 



a. 



e'ij , 
3Ke, 



(63) 



where the real parameters fi and K are the shear and bulk moduli, respectively. Devia- 
toric stresses and strains (physical components) are defined by 







1 ~A 


^ij 


= CTij 




~l 






^ij 


= eij 


- 3 e (>ij > 



(64) 



where a and e are the traces of the stress and strain tensors. 



(T 

e 



(65) 



with implicit summation over repeated indices. 

The correspondence principle \Fung 1965 states that the constitutive equations for a 
viscoelastic body are similar in form but with complex shear modulus fl{oo) and complex 
bulk modulus K{ijj) depending on the angular frequency uj: 



^'ij 



a 



^ij 



2/i(w) 
3K{uj)e 



(66) 
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The dependence on frequency of and K(uj) defines the rheology of the material. 



Peale and Cassen |1978j gave the correct formula for the rate at which stresses do 
work per unit volume: 

P{t) = a^J{t) e^J{t) . (67) 



Kaula\ [1964 incorrectly ntroduced a factor 1/2 in this formula (see his Eqs. (10)-(16)- 



(19)), arguing that 'work is done only with motion, not with just change of force'. 
Moreover Kaula thought that a fraction of this energy (1/2 for the Moon, 1/27 for the 
Earth) went into orbital kinetic energy. 

The dissipation power per unit volume averaged over one orbital period T reads 



P 



1 

f 



aij{t) eij{t) dt . 



(68) 



In terms of the Fourier components appearing in Eqs. (62), the dissipated power reads 

(69) 



P = ~Im(a^je*j) , 



where the asterisk denotes the complex conjugate. If stresses and strains are in phase, 
Ojj^^j^Jsjceal so that the dissipated power vanishes. Eq. (69) agrees with Eq. (24) of 



Tobie et al. 2005b . 



In terms of deviatoric tensors, the dissipated power reads 

1 



P 



}Im{jl) e'i* + -Im{K) |e| 



(70) 



The first term is due to the change of body shape whereas the second term is due to 
the change of volume. Note that Eq. (10) of Kaula [1964 has an incorrect factor 1/6 
instead of 1/2 in the term proportional to K. 

In terms of non-deviatoric strains, the dissipated power reads 



P = LoIm{ji) 



(71) 



Note that Eq. (8) oi \Segatz et aL\ [1988 is wrong by having a factor uj/2 instead of uj in 
front of Im[fi). 

'Fo\\owmg \Ojakangas and Stevenson 1989a , many authors write the dissipated power 
in terms of the averaged 'squared' strain rate defined as 



1 
T 



hj{t) €ij{t) dt 
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For an incompressible body with Maxwell rheology, Eq. (71) becomes 



P 



;2 



(73) 



oj 1 + (cotm)^ 

where tm = is the Maxwell time {rj is the viscosity). Recall that Maxwell rheology 
is defined by [e.g. \Peltier\ |1982[ \Tobie et al.\ |2005bt | Wahr et al.\ |2009| 



{UJTM + i) 



1 + (wrM)2 

Appendix B. Differential operators on the sphere 



(74) 



The tidal potential at the surface of a spherical planet is a scalar function which means 
that it is invariant under rotations of the coordinate system. The dissipated power 
at a given radius is also a scalar function which depends on derivatives of the Fourier- 
transformed tidal potential. My purpose is to construct a scalar function from derivatives 
of another scalar function. Let f(6, (p) be a complex scalar function defined on the sphere. 
How can we combine derivatives of / into a scalar? The standard method consists in (1) 
constructing tensors with covariant derivatives Vq and (2) contracting all indices with 
the inverse metric on the unit sphere, 



g'^P = diag{l, sin 



(75) 



with Greek indices denoting angular coordinates {O^cj)). Covariant derivatives and other 
operators on the sphere are discussed in detail in Appendices B and D of Beuthe [2008| . 

A well-known scalar expression of order two in derivatives is the contracted double 
covariant derivative, called the spherical Laplacian: 



A/ = 5"^V«V./. 



(76) 



The eigenfunctions of A are spherical surface harmonics of degree i with eigenvalues 

5e = -£{£ + 1) . (77) 
Contracting covariant derivatives of / and /* yields scalars of order two and four: 



^?2(/,r) 



g'^'g^' (v,v^/) (v.v^r 



In terms of usual derivatives, Vo reads 



didr 

de de ' sin^e 



+ 



1 df df 



(78) 



(79) 
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Since P4 has a complicated expression in terms of usual derivatives, I define auxiliary 
operators as follows. In the normalized basis Beuthe 2008 , Section 3.2 and Appendix A], 
the operators corresponding to double covariant derivatives are given by 



Oif 



— 2^ V^V^/, 



sm 
1 



sin( 



(80) 



The operators Oi are given in terms of usual derivatives by 

52 



Oi 
O3 



89'^ ' 
1 

sin^ 6 
1 

sin 9 



\d9dip 



+ cos 9 sin ( 



do 



cot 9 



d_ 

dip 



(81) 



The operators A and P4 can be expressed in terms of usual derivatives through combi- 
nations of the operators Oi: 



A/ 



(01 + 02)/, 

(Oi/) (OiD + (O2/) (02/*) + 2 (O3/) (OsD • (82) 



Finally, the operators 2?2 and P4 expressed in terms of covariant derivatives (Eqs. (78)) 
can be transformed with Eqs. (G.2)-(G.3) of Beuthe [2010 so that derivatives appear 
only through the Laplacian. If / is a spherical harmonic of degree i, I get 



'D2if,n 



'A-25,)\f\\ 

AA - 2 {26i + 1) A + 46i {6i + l))\f\ 



(83) 



where 5e results from the action of A on / or /* (see Eq. (77)). However |/| is not 
an eigenfunction of A because it is the sum of harmonic components of even degrees 
between and 2i. 



Appendix C. Tidal potential for synchronous orbit 



Let {R, 9, (j)) be the radius, colatitude and longitude of point on the surface of a satellite 
orbiting with semi-major axis a around a point mass M (the primary). The time- 
dependent part of the tidal potential at {R, 9,(p) is given at first order in eccentricity e 
and obliquity / by Eq. (1) of Kaula 19641: 



U{t,9,(t>) 



GMR^ ( 3e 



P2o(cos6') COS 
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+ o A2(cos(9) (- COS (A^ - 2A + 2ujp) + 7cos {3M - 2A + 2ujp)) 
8 

+ ^sin/P2i(cos6')(sm(2X-A + 2wp) + sinA)^ . (84) 

Pim{x) are the associated Legendre functions of degree £ and order m, Ai is the mean 
anomaly and cOp is the argument of pericentre. The auxiliary variable A is defined by 

A = (j) + e-Q, (85) 

where (j) is the longitude in a frame (x, y, z) attached to the satellite (x and y are within 
the equatorial plane) , Q is the sidereal time for the reference meridian (the angle between 
the X-axis and an inertially fixed point Q on the equator) and O is the longitude of the 
ascending node (measured from Q). 

If the orbital motion is synchronous with the rotation of the satellite, the sidereal 
time and the mean anomaly are related by = + const. Note that Ai gives the 
position of the planet only at pericentre and apocentre, unless eccentricity vanishes, and 
that the angles {Q,Ai) are not measured in the same plane, unless obliquity vanishes. I 
choose the direction of the x-axis by setting 

Q = M+vj. (86) 

where w is the longitude of pericentre: 

■CD = ^} + Up . (87) 

If the obliquity vanishes, the x-axis points to the primary when the latter is at pericentre: 
Q = zu when = {mod 27r). If the eccentricity vanishes, the x-axis points to the 
primary when the latter is at the ascending node: Q = 0, when M +ujp = {mod 2tt). In 
the general case of non-zero obliquity and eccentricity, the x-axis points to the primary 
neither at pericentre nor at the ascending node. 

Substituting Eq. (86) into the tidal potential (84) and setting M. = nt (time of 
pericentre passage set to zero, n being the mean motion), I get 

,2/ 3 „ ^ , /3 



U {t, 9, (j)) = {nR) y~ 2^ -f2o('^os ^) cosnt + e P22(cos ^) ^ "^^^ ""^^ 2(j) + sin nt sin 20 

+ sin I P21 (cos 9) sin {nt + uop) cos < 

The approximation GM/a? = is valid if the deformed body is much lighter than the 
primary. 

The Fourier coefficient of the tidal potential defined by Eq. ^ is thus equal to 

= (nP)^(^-^eP2o(cos6') + eP22(cos6') cos2(/) - z sin 20^ - i e*'^*' sinIP2i(cos6') cos( 

(89) 

The part of the tidal potential due to non-zero obliquity is dephased by the argument 
of pericentre Up (frame-independent quantity) with respect to the part due to non-zero 
eccentricity. 
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Appendix D. Squared norm of the tidal potential 

Let the Fourier coefficient of an arbitrary tidal potential of degree two be written as 

= {nRf (a P2o(cos 6) + (/3 cos 2</> + 7 sin 20) P22(cos6') + (6 cos0 + esin0) P2i(cos^)) , 

(90) 

where (a, /3, 7, 5, e) are dimensionless complex numbers. The term in P20 appears when 
the orbit is eccentric, the terms in P22 appear either when the rotation is not synchronous 
with the orbital motion or when the orbit is eccentric or when there are forced librations 
in longitude, whereas the term depending on P21 appears when the obliquity is not zero. 

The squared norm of the tidal potential involves products of Legendre functions of 
degree two. These can be expressed as linear combinations of Legendre functions with 



the formulas of Balmino 1994, pp. 358-359] to which I add 



P21P22 = (3/35) P43 = (12/7) P21 - (18/35) P41 • (91) 



The coefficients {(lim^^im) even order in the spherical harmonic expansion (34) 
read 



'^00 


= (l/5)(|ap + 12(|/3|2 + |7|2)+3(|<5|2 




0-20 


= (l/7)(2|a|2-24(|/3|2 + |7|2)+3(|5| 




0-22 


= (1/14) {-8Re[al3*] + 3 {\5\^ - \e\^)) , 




^22 


= (l/7)Pe[-4a7* +3(5e*] , 




a^Q 


= (l8/35)(|a|2 + 2(|/3|2 + |7|2-|5|2- 




042 


= (3/35) {2Re[a(3*] + |(5|2 - jep) , 




h2 


= {6/35) Re[a-f* + 6 e*], 




0-44 


= (3/70)(|/3|2-|7|2) , 




644 


= (3/35) Pe[/37*] , 





(92) 

while the coefficients of odd order (resulting from interferences between non-zero obliq- 
uity terms and other terms) are given by 

= (2/7) Pe[(a + 6/3)5* + 676*], 

= (2/7) Re[{a - 6/3)e* + 6-f6*] , 

041 = (18/35) Pe[(a - /3)(5* - 76*] , 

641 = (18/35) Pe[(a + /3)e* - 7(5*] , 

= (3/35)Pe[/3(5* -7e*], 

= (3/35) Pe[/3e* + 7(5*] . (93) 

Appendix E. Variance of spatial patterns 



The variance of a function ^ defined on the unit sphere reads 

1 



var («') = -'-/ (^' - ^')^dcj, (94) 
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where ^ is the average of \I' on the sphere and da is the integration measure. The 
angular functions ^ (j take the generic form 

^j = x^Q + y^^ + z^^ {J = A,B,C). (95) 



with coefficients {x,y,z) given by Eqs. (36). The variance of ^'j is equal to 



var(^'j) = var(4'2) + z'^ var(^'4) , (96) 
where each term can be evaluated with the expansion (|34|), 



var (M/,) = Vim + (beJ^) ■ (97) 

m=0 

The coefficients Vim are computed with the orthogonality relation of spherical harmonics: 

1 3 12^ 



(^20,^21,^^22) = 1 g 1 , 

(y40,^41,^42,'^43,'^44) = ^ , y , 20, 280, 2240^ . (98) 

Table |3| gives the standard deviation (square root of the variance) of ^ ^ for syn- 
chronous rotation (eccentricity tides only or obliquity tides only). 



Appendix F. 
density 



Multi-layer incompressible body of uniform 



Consider an incompressible body of radius R stratified into spherical homogeneous layers 
which can be solid or liquid, except the surface layer which is solid. The yi functions 
(defined as in Takeuchi and Saito 1972 ) required for the strains (Eqs. (|8|-([9])) have the 



following form in the solid layers Sabadini and Vermeersen , 2004 : 



h2 

9 
h2 
9 

9 



5 -3 1 . - 

— ar -\ — or 

42 2 



dr 



21 



ar^ + br + - cr + 



dr 



(99) 



where f = r/R, g is the surface gravity, fx is the viscoelastic shear modulus of the layer 
and /i2 = gyi{R) is the tidal Love number for radial displacement. The dimensionless 
coefficients {a,b,c,d) are constant in each layer. The factorization of /12/5' introduces 
one more unknown (/i2) in the problem but also one more constraint: 



1,1, 

-as + bs + -cs + ds 



1 



(100) 
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where the index S denotes the surface layer. The boundary condition of no shear stress 
at the surface {U/^iR) = 0) yields another constraint: 



— as + bs + -cs + -ds = 0. 



(101) 



The two other boundary conditions of the problem involve the gravitational potential so 
that it is necessary to solve for strains, stresses and gravitational potential at the same 
time. 

Suppose now that the density p of the body is uniform. In that case, it is possible 
to solve for the deformations independently of the gravitational potential. First, it can 
be shown that the dimensionless constants (a, b, c, d) depend only on the ratios between 
the parameters characterizing the different layers, that is their radii and shear moduli 
(the proof relies on the fact that the transition matrix between the layers i and i + 1 
depends only on Ri/R and fli/fli+i). As a result the functions {yi,y^,ry^/ jl) depend on 
a dimensional viscoelastic parameter only through the prefactor h2 ■ Therefore the ratios 
of squared strains (determining the relative importance of the dissipation patterns) do 
not depend on the absolute scale of deformation (parameterized by jls/{pgR)) but only 
on the ratios between the parameters characterizing the different layers (radii and shear 
moduli). 

Second, the boundary condition on the gravitational potential given by Ry^^R) = 5 
\Takeuchi et al. 1962; Saito, 1974] ) implies that 



ho 



(102) 



where the tidal Love number for gravity is defined by k2 = y^{R) — 1 (2/5 is the gravita- 
tional potential) . This relation was expected since the radial deformation of the surface 
completely characterizes the gravity potential of the body when the density is uniform. 
Third, the boundary condition on the radial stress {y^iR) = 0) combined with 



Eqs. (lOO)-(lOl) implies that 



35 \ PS 



-1 



Fourth, the tidal Love number I2 = gy^{R) is related to /12 by 



3 

10 



cs 



6 



ds]h 



(103) 



(104) 



In Eqs. (103)-(104), the constants {cs,ds) must be determined by relating them with 



propagation matrices to two unknown constants in the core (if it is solid) and applying 
the conditions (lOO)-(lOl). Each liquid layer introduces two additional unknowns but 

[e.g. 



also two additional constraints 



Vermeersen . 2011 



Sabadini and Vermeersen, 2004; Jara-Orue and 
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For a one-layer body, the constants {c,d) in Eqs. (99) vanish otherwise the solution 

3ly yiel 

19 /i 



diverges in r = 0. Eqs. (103)-(104) immediately yield the well-known formulas 

-1 



ho 



1 + 



2 pgR 



In that case, the solution of Eqs. ( 100)-( 10l| ) is 

{a,b,c,d) = (-21/5,8/5,0,0) , 
yielding the displacement functions 



2/1,2/3, 



ho r 



2. 5^2^g^^__2^ 



(105) 



(106) 



(107) 



Appendix G. Two-layer incompressible body of uniform den- 
sity 

Consider an incompressible body of uniform density p composed of a viscoelastic core 
of radius Rc and a viscoelastic mantle of radius R with shear moduli pc and pM, 
respectively. The surface gravity is denoted g. I define the dimensionless ratios 



(r , X , ^ , Pred) 



Rc fj-c fj-M 



R R PM pgR 



(108) 



The coefficients (o, b, c, d) in the displacement functions ( |99[ ) are ratios of polynomials 
in X (dimensionless core radius) and (ratio of shear moduli): 



{a,b,c,d) = {pa,Pb,Pc,Pd) /q- 
The polynomial in the denominator is given by 

q = 4(24 + 40x3-45x'^-19x^°) + 2(89 + 15x3-5x^ + 76x^°)C 
+ 38 (2 - 5 + 5 x"^ - 2 x^°) . 

In the core, Pc = Pd = while pa and pb are given by 



(109) 



(110) 



Pa 
P? 



-294 



5+ (2-5x3 + 8x^) (^-1) 



280 + (152 - 21 x^ - 147 x^ + 296 x^) (^ - 1) 



(111) 



In the mantle, pc and pd are given by 



P^ 



Pf 



x^ {i - 1) [l6 (40 - 21 x^ - 19 x'^) + 19 (40 - 21 x^ + 16 x'^) i 



x^ (^ - 1) [2 (64 - 45 x^ - 19 x^) + 19 (8 - 5 x^ + 2 x^) ^ 



:il2) 
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while Pa and are related to them by Eqs. (lOO)-(lOl): 



-(21/5) (?-7p*^ 
(8/5)(?-(l/2)p, 



M 
c 



(113) 



The substitution of Eqs. (112) into Eqs. (103)-( [T04 ) yields the displacement tidal Love 
numbers (^2 is given by Eq. (102)): 



h2 
h 



1 + 



19 

y 



Ph 



-1 



3 

10 



(114) 



with 



Ph 
Pi 



lOpf + (35/2)pf , 
(l/4)pf + (7/6)pf . 



(115) 



If 1^1 = (liquid core) or if |^| — t- oo (infinitely rigid core), viscoelastic parameters 
appear in the functions yi through the prefactor /i2 so that tidal deformations depend 
depend only on the core radius x. The Love number /12 becomes 



/l2(x) = 

If the core is liquid, z{x) is given by 



^(^1 + z{x) fJred 



z"''(x) = 12 



19 - 75 x3 + 112 - 75 x^ + 19x 

19x10 



10 



24 + 40 x3 - 45 x7 
If the core is infinitely rigid, z(x) is given by 

38 + 225x3 - 336x5 _^ 200x^ + 48 x^^ 



z"f(x) 



2(2-5x3 + 5x^-2x10) 



(116) 



(117) 



(118) 



The thin shell limit is obtained first by setting ^ = (liquid core) and then by taking 
the limit x — >• 1. In that case, z'*'^(x) ~ (60/ll)(l — x) and pi/q — t- 3/110 so that 



h2 



5 
2 
3 

n 



-, 60 ,^ 



l'2 ■ 



(119) 
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